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Abstract: Contributions of the KASCADE-Grande Collaboration to the 32nd International Cosmic Ray Conference, 
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Abstract: The detection of high-energy cosmic rays above a few hundred TeV is realized by the observation of extensive 
air-showers. By using the multi-detector setup of KASCADE-Grande, energy spectrum, elemental composition, and 
anisotropies of high-energy cosmic rays in the energy range from below the knee up to 1 EeV are investigated. The most 
distinct feature of the spectrum, the knee, is thought to be the beginning of the end of the galactic origin of cosmic rays. 
As the highest energies (above the ankle) are most probably of extragalactic origin, between 10 PeV to 1 EeV one expects 
the transition of galactic to extragalactic origin. KASCADE-Grande is dedicated to explore this transition region. The 
estimation of energy and mass of the high-energy primary particles is based on the combined investigation of the charged 
particle, the electron, and the muon components measured by the detector arrays of Grande and KASCADE. The latest 
analysis results have shown that a knee-like structure in the all-particle energy spectrum at 83 PeV is due to a decrease 
of flux of the heavy mass component. In this contribution an overview is given on various different analysis methods to 
verify this finding. 

Keywords: KASCADE-Grande, 10-1000PeV, spectrum and composition 



KASCADE-Grande: Main parts of the experiment are the 
Grande array spread over an area of 700 x 700 m 2 , the 
original KASCADE array covering 200 x 200 m 2 with un- 
shielded and shielded detectors, and additional muon track- 
ing devices. This multi-detector system allows us to inves- 
tigate the energy spectrum, composition, and anisotropies 
of cosmic rays in the energy range up to 1 EeV. The estima- 
tion of energy and mass of the primary particles is based on 
the combined investigation of the charged particle, the elec- 



tron, and the muon components measured by the detector 
arrays of Grande and KASCADE. 

The multi-detector experiment KASCADE [1] (located at 
49.1°n, 8.4°e, HOma.s.l.) was extended to KASCADE- 
Grande in 2003 by installing a large array of 37 stations 
consisting of 10 m 2 scintillation detectors each (fig. 1). 
KASCADE-Grande [2] provides an area of 0.5 km 2 and 
operates jointly with the existing KASCADE detectors. 
The joint measurements with the KASCADE muon track- 
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Figure 1: Layout of the KASCADE-Grande experiment: 
The original KASCADE (with Array, Muon Tracking De- 
tector and Central Detector), the distribution of the 37 sta- 
tions of the Grande array, and the small Piccolo cluster for 
fast trigger purposes are shown. The outer 12 clusters of 
the KASCADE array consist of [i- and e/7-detectors, the 
inner 4 clusters of e/7-detectors, only. 

ing devices are ensured by an additional cluster (Piccolo) 
close to the center of KASCADE-Grande for fast trigger 
purposes. For results of the muon tracking devices see ref- 
erences [3, 4]. While the Grande detectors are sensitive to 
charged particles, the KASCADE array detectors measure 
the electromagnetic component and the muonic component 
separately. These muon detectors enable to reconstruct the 
total number of muons on an event-by-event basis also for 
Grande triggered events. 

Basic shower observables like the core position, angle-of- 
incidence, and total number of charged particles are pro- 
vided by the measurements of the Grande stations. A 
core position resolution of « 5 m, a direction resolution of 
w 0.7°, and a resolution of the total particle number in the 
showers of sa 15% is achieved. The total number of muons 
(N/x resolution « 25%) is calculated using the core posi- 
tion determined by the Grande array and the muon densities 
measured by the KASCADE muon array detectors. Full ef- 
ficiency for triggering and reconstruction of air-showers is 
reached at primary energy of « 10 16 eV, slightly varying 
on the cuts needed for the reconstruction of the different 
observables [2]. 

The strategy of the KASCADE-Grande data analysis to 
reconstruct the energy spectrum and elemental composi- 
tion of cosmic rays is to use the multi-detector set-up of 
the experiment and to apply different analysis methods to 
the same data sample. This has advantages in various as- 
pects: One would expect the same results by all methods 
when the measurements are accurate enough, when the re- 
constructions work without failures, and when the Monte- 




3.5 4 4.5 5 5.5 6 6.5 7 7.5 

muon number log (N ) 



Figure 2: Two-dimensional distribution of the shower sizes 
charged particle number and total muon number as mea- 
sured by KASCADE-Grande. All quality cuts are applied, 
i.e. these data are the basis for the mass composition stud- 
ies. 

Carlo simulations describe correctly and consistently the 
shower development and detector response. 

The main air-shower observables of KASCADE-Grande, 
shower size and total number of muons, could be recon- 
structed with high precision and low systematic uncertain- 
ties and are used in the following for the data analysis 
(fig. 2). 

The all-particle energy spectrum: In a first step of the 
analysis, we reconstructed the all-particle energy spec- 
trum. Applying various reconstruction methods to the 
KASCADE-Grande data the obtained all-particle energy 
spectra are compared for cross-checks of the reconstruc- 
tion, for studies of systematic uncertainties and for testing 
the validity of the underlying hadronic interaction models. 
By combining both observables and using the hadronic in- 
teraction model QGSJet-II, a composition independent all- 
particle energy spectrum of cosmic rays is reconstructed in 
the energy range of 10 16 eV to 10 18 eV within a total un- 
certainty in flux of 10-15%. 

Despite the overall smooth power law behavior of the re- 
sulting all-particle spectrum, there are some structures ob- 
served, which do not allow to describe the spectrum with 
a single slope index [5]. There is a feature in the spec- 
trum showing a small break at around 8 • 10 16 eV. The 
power law index of 7 = —2.95 ± 0.05 is obtained by fit- 
ting the range before this break. Applying a second power 
law above the break an index of 7 — —3.24 ± 0.08 is ob- 
tained (see figure 3). With a statistical significance of more 
than 2 sigma the two power laws are incompatible with 
each other. This slight slope change occurs at an energy 
where the rigidity dependent knee of the iron component 
would be expected (KASCADE QGSJet based analysis as- 
signs the proton knee to an energy of 2 — 4 • 10 15 eV). De- 
spite the fact, that the discussed spectrum is based on the 
QGSJet-II hadronic interaction model, there is confidence 
that the found structures of the energy spectrum remain sta- 
ble. Tests with EPOS as well as investigations of the spec- 
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Figure 3: Reconstructed all -particle energy spectrum to- 
gether with the spectra of the electron-poor and electron- 
rich components. Fits on the spectra and resulting slopes 
are also indicated. 

tra of the pure (and independently obtained) observables 
shower size and total muon number, confirmed the struc- 
tures [5]. 

Composition: A conclusion on the origin of the found 
structures in the all-particle spectrum is not possible with- 
out investigating the composition in detail in this energy 
range. The basic goal of the KASCADE-Grande experi- 
ment is the determination of the chemical composition in 
the primary energy range 10 16 — 10 18 eV. Like for the re- 
construction of the energy, again several methods using 
different observables are applied to the registered data in 
order to study systematic uncertainties. However, the in- 
fluence of predictions of the hadronic interaction models 
has a much larger influence on the composition than on 
the primary energy. As it is well known from KASCADE 
data analysis [6] that the relative abundances of the indi- 
vidual elements or elemental groups are very dependent on 
the hadronic interaction model underlying the analyses the 
strategy is to derive the energy spectra of the individual 
mass groups. The structure or characteristics of these spec- 
tra are found to be much less affected by the differences 
of the various hadronic interaction models than the relative 
abundance. The present goal is to verify the structure found 
in the all-particle energy spectrum at around 100 PeV in the 
individual mass group spectra and to assign it to a particu- 
lar mass. 

The main observables taken into account for composition 
studies at KASCADE-Grande are the shower size (N c h, 
or the subsequently derived electron number N e ) and the 
muon shower size (N^). Figure 2 displays the correlation 
of these two observables, i.e. this distribution is the basis 
of the composition analysis with KASCADE-Grande data. 
For all the methods it is crucial to verify the sensitivity of 
the observables to different primary particles and the repro- 
ducibility of the measurements with the hadronic interac- 
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Figure 4: Shower size ratio distributions for a certain bin 
in charged particle number. Shown is the measured distri- 
bution as well as the simulated ones for three primary mass 
groups and the resulting sum. 

tion model in use as a function of sizes and the atmospheric 
depth. So far, in the composition analysis we concentrate 
on interpreting the data with the hadronic interaction model 
QGSJet-II (and FLUKA as low energy interaction model). 

Four methods of composition studies at KASCADE- 
Grande are discussed in the following, all showing the 
same result: The knee-like structure in the all-particle spec- 
trum is due to the decrease of the flux of the heavy compo- 
nent of primary cosmic particles. 

• Charged particle - muon number ratio [7]: The total 
number of charged particles N c h and the total num- 
ber of muons of each recorded event are con- 
sidered and the distribution of N^/N c h is studied in 
different intervals of N c h (corresponding to different 
energy intervals) and zenith angle. The experimen- 
tal distribution of the observable N^/Nch is taken 
into account and fitted with a linear combination of 
elemental contributions from simulations, where we 
distinguish three groups: light, medium and heavy 
primaries. By this way the means and the widths of 
the distributions as two mass sensitive observables 
are taken into account (see as an example figure 4). 
The width of the data distribution is in all ranges 
of N c h so large that always three mass groups are 
needed to describe them. Using the Monte Carlo 
simulations the corresponding energy of each mass 
group and shower size bin is assigned to obtain the 
spectra of the individual mass groups. Investigating 
the spectra of these individual mass groups, it is ob- 
vious that the break seen in the all particle spectrum 
is due to heavy primary masses. 

• The Y-cut method [8]: Here, the shower ratio 
Ycic = logiV^/ logiVcTi between the muon and 
the charged particle numbers, both corrected for at- 
mospheric attenuation by the Constant Intensity Cut 
method (CIC), is used as parameter to separate the 
KASCADE-Grande data into different mass groups. 
MC simulations performed with CORSIKA on the 
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Figure 5: KASCADE and KASCADE-Grande recon- 
structed energy spectra of individual mass groups. 

framework of FLUKA/QGSJET-II are employed to 
obtain the expected Ycic distributions as a function 
of the energy for different cosmic ray primaries as a 
basis for the separation. Then the Ycic -parameter 
is used to divide the KASCADE-Grande data into 
electron-rich and electron-poor events, i.e. generated 
by light and heavy primaries. The results confirm the 
findings above. 

• The k-parameter method [9]: Using the above men- 
tioned reconstruction of the energy spectrum by cor- 
relating the size of the charged particles N c h and 
muons on an event-by-event basis, the mass 
sensitivity is minimized by means of the parameter 
k(N c h, Nfj). On the other hand, the evolution of k as 
a function of energy keeps track of the evolution of 
the composition, and allows an event-by-event sepa- 
ration between light and heavy primaries. Using k 
as separation parameter for different mass groups, 
where the exact values of k have to be determined 
with help of simulations, directly the energy spectra 
of the mass groups are obtained. Figure 3 shows the 
resulting spectra. As the significance of the break of 
the all-particle spectrum increases significantly when 
enhancing heavy primaries, it is obvious that this 
structure is introduced by heavy primaries. 

• The Gold-unfolding method [10]: This method 
is based on the unfolding of the two-dimensional 
shower size spectrum (fig. 2) in a similar way as it 
was developed for the KASCADE data analysis [6]. 
Due to the fact that the accuracies in reconstructing 
the shower sizes for KASCADE-Grande are not as 
high as in case of KASCADE, for comparing both 
results, primary masses are combined to light (H), 
medium (He+C+Si) and heavy (Fe) spectra. The re- 
sulting individual mass group spectra can be com- 
bined to provide a solution for the entire energy 



range from 1 PeV to 1 EeV. This analysis was per- 
formed for the combination of N c h with [10] as 
well as for the combination N e with N^. Figure 5 
shows the results for the latter case [11]. The ob- 
tained spectra are compared with the results using 
KASCADE data in the lower energy range, both us- 
ing the hadronic interaction model QGSJet II with 
FLUKA. The agreement in the overlapping energy 
range of the spectra is remarkable. The knee-like 
feature in the iron component is found at an energy 
which is about 26 times higher than the proton knee 
in case of KASCADE. 

Summarizing, by these first composition studies it is seen 
that the knee-like feature in the energy spectrum at about 
8 • 10 16 eV is due to a kink in the spectrum of the heavy 
component of primary cosmic particles. In addition, at 
least three primary elements are needed to describe the ex- 
perimental data over the entire energy range accessible by 
KASCADE-Grande, i.e. up to 1 EeV. Furtheron, it was 
found that QGSJet-II, the hadronic interaction model in 
use, can fairly well reproduce the data and, in particular, 
provides a consistent solution on the elemental composi- 
tion, independent of the method in use. One has to remark, 
that using another hadronic interaction model would proba- 
bly lead to significant changes in the relative abundances of 
the elemental groups as different models predict different 
shower sizes for a certain energy and mass of the primary 
cosmic ray. But, we are confident that the obtained spectral 
form for the heavy and light component of the cosmic ray 
spectrum will remain unchanged. 
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A study of the mass composition of cosmic rays based on an event-by-event assignment with 
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Abstract: The cosmic ray energy spectrum between 10 16 eV and 10 18 eV is reconstructed in KASCADE-Grande by 
correlating the size of the charged particles (N c /j) and muons (N M ) on an event-by-event basis. In the energy assignment, 
the mass sensitivity is minimized by means of a parameter k(N c h, N^). On the other hand, the evolution of k as a 
function of energy keeps track of the evolution of the composition, and allows an event-by-event separation between 
electron rich and electron poor primaries. A first result on the evolution of k and its connection with the shape of the 
energy spectrum, is presented in the framework of the QGSJet 11-03 interaction model. 

Keywords: Composition, energy spectrum, 10 16 - 10 18 eV, KASCADE-Grande. 



1 Introduction 

The study of the energy spectrum and of the chemical com- 
position of cosmic rays are fundamental tools to understand 
origin, acceleration and propagation of cosmic rays. The 
energy range between 10 16 eV and 10 18 eV is quite impor- 
tant from astrophysical point of view because it is expected 
that in this energy range the transition between galactic 
and extra-galactic origin of cosmic rays will occur. The 
results obtained at lower energies by KASCADE [1] and 
EAS-TOP [2] as well as by other experiments suggest that 
the knee in the primary energy spectrum around 3 - 4 x 



10 15 eV is due to the break in the spectra of light elements. 
Moreover, KASCADE results seem to indicate that a rigid- 
ity dependent mechanism is responsible for such knees. 
Therefore, a knee of the heaviest components would be ex- 
pected in the range of 10 16 eV to 10 18 eV. Various theories 
with different assumptions try to explain the rather smooth 
behavior of the cosmic ray energy spectrum in this energy 
range (i.e. [3, 4]). In order to discriminate between the dif- 
ferent models, a very precise measurement of the possible 
structures of the energy spectrum and of the evolution of 
the composition is needed. 
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2 The Technique 

The technique employed to derive the all-particle energy 
spectrum and the abundance of electron-rich (e.r.) and 
electron-poor (e.p.) primaries is based on the correlation 
between the size of the charged particles (N c /j) and muons 
(N M ) on an event-by-event basis. The method itself has 
been described in detail in [5]. 

A sample of Monte Carlo data was simulated including 
the full air shower development in the atmosphere, the re- 
sponse of the detector and its electronics as well as their un- 
certainties. In this way, the parameters reconstructed from 
simulation are obtained in the same way as for real data. 
The EAS events were generated with an isotropic distribu- 
tion with spectral index 7 = — 3 and were simulated with 
CORSIKA [6] and the hadronic Monte Carlo generators 
FLUKA [7] and QGSJet 11-03 [8]. Sets of simulated events 
were produced in the energy range from 10 15 eV to 10 18 eV 
with high statistics and for five elements: H, He, C, Si and 
Fe, representative for different mass groups (« 353.000 
events per primary). Few events up to 3 • 10 18 eV were also 
generated in order to cross-check the reconstruction behav- 
ior at the highest energies. 

Grande stations [9] are used to provide core position and 
angle-of-incidence, as well as the total number of charged 
particles in the shower, by means of a maximum likeli- 
hood procedure comparing the measured number of par- 
ticles with the one expected from a modified NKG lateral 
distribution function [10] of charged particles in the EAS. 
The total number of muons is calculated using the core po- 
sition determined by the Grande array and the muon den- 
sities measured by the KASCADE muon array detectors. 
The total number of muons in the shower disk (above 
the energy threshold of 230 MeV) is derived from a max- 
imum likelihood estimation assuming a lateral distribution 
function based on the one proposed by Lagutin and Raikin 
[11]. The reconstruction procedures and obtained accu- 
racies of KASCADE-Grande observables are described in 
detail in reference [9]. 

For the reconstructed events, we restricted ourselves to 
events with zenith angles lower than 40°. Additionally, 
only air showers with cores located in a central area on 
KASCADE-Grande were selected. With this cut on the 
fiducial area, border effects are discarded and possible 
under- and overestimations of the muon number for events 
close to and far away from the center of the KASCADE 
array are reduced. All of these cuts were applied also to 
the Monte Carlo simulations to study the effects and to op- 
timize the cuts. Full efficiency for triggering and recon- 
struction of air-showers is reached at primary energy of 
py 10 16 eV, slightly depending on the cuts needed for the 
reconstruction of the different observables [9]. 
The analysis presented here is finally based on 1173 days 
of data and the cuts on the sensitive central area and zenith 
angle correspond to a total acceptance of A = 1.976 • 10 9 
cm 2 - sr, and an exposure of N = 2.003 • 10 17 cm 2 - s-sr, 
respectively. 

With Monte Carlo simulations a formula is obtained to 
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Figure 1: Scatter plot of the reconstructed N^/N^ vs. N c h 
for H and Fe primaries for the first angular bin. The full 
dots and error bars indicate the mean and its statistical error 
of the distribution of the individual events (small dots). The 
fits result in parameters c and d of expression 3. 




Figure 2: Scatter plot of E vs. N c h for Fe and H primaries. 
The fits result in parameters a and b of expression 1 . 

calculate the primary energy per individual shower on the 
basis of N c h an d N^. The formula takes into account the 
mass sensitivity in order to minimize the composition de- 
pendence in the energy assignment, and at the same time, 
provides an event-by-event separation between e.r. and e.p. 
candidates. The formula is defined for 5 different zenith an- 
gle intervals (9 < 16.7, 16.7 < 9 < 24.0, 24.0 < 9 < 29.9, 
29.9 < 9 < 35.1, 35.1 < 9 < 40.0) independently, to take 
into account the shower attenuation in atmosphere. Data 
are combined only at the very last stage to obtain a unique 
power law spectrum and mass composition. The energy 
assignment is defined as E = f{N c h, k) (see equation 1), 
where N c /j is the size of the charged particle component 
and the parameter k is defined through the ratio of the sizes 
of the N c h and muon (N M ) components: k = g(N c f l ,N fl ) 
(see equation 2). The main aim of the k variable is to take 
into account the average differences in the N C ^/N M ratio 
among different primaries with same N c jj and the shower to 
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Table 1: Parameters of the calibration functions. 



Angles [deg] 


a 


b 


c 


d 




H 


Fe 


H 


Fe 


H 


Fe 


H 


Fe 


0.0 - 16.7 


0.910 


0.876 


1.333 


1.817 


0.100 


0.161 


0.786 


-0.055 


16.7 - 24.0 


0.894 


0.878 


1.495 


1.923 


0.081 


0.179 


0.884 


-0.254 


24.0 - 29.9 


0.937 


0.889 


1.301 


1.935 


0.104 


0.156 


0.677 


-0.170 


29.9 - 35.1 


0.934 


0.881 


1.458 


2.099 


0.109 


0.171 


0.543 


-0.351 


35.1 - 40.0 


0.919 


0.875 


1.748 


2.287 


0.105 


0.156 


0.412 


-0.348 



shower fluctuations for events of the same primary mass: 



log 10 {E[GeV}) = 

logxo(N ch /N„) - log w (N ch /N^ H 



k = 



an + (a>Fe - a H ) ■ k] ■ log w (N ch ) 
h b H + (bFe -b H )-k 



(2) 
(3) 



log w (N ch /N„) Fe - logwiNch/N^H 
logwiNch/N^H.pe = c H ^ Fe ■ log 10 {N ch ) + d H , Fe 
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Figure 3: Ratio between the reconstructed and true simu- 
lated energy spectra for H, Fe and all mixed primaries sum- 
ming up all angular bins. 

defined for H and Fe in the first angular bin, while Fig. 3 
shows the capability of reproducing simulated energy spec- 
tra. Pure spectra of H, Fe and a mixture of 5 different pri- 
maries with 20% abundance each are shown as examples. 
The true flux is always reproduced within 10% uncertainty. 

3 Results 

Fig. 4 shows, for the sum of the first two angular bins, the 
evolution of the k parameter as a function of the recon- 
structed energy obtained by simulating with QGSJet 11-03 
model pure H, He, C, Si, Fe primary spectra, as well as for 
the experimental data. Similar behavior is obtained also for 
the other angular bins. The error bars indicate the average 



(1) 



The k parameter is, by definition of eq. (2), a number cen- 
tered around for H showers and 1 for Fe ones if expressed 
as a function of N c h, while slightly shifted when reported 
as a function of energy (see Fig. 4). The complete list of 
parameters a -d can be found in table 1 . 

Figs. 1, 2 show the scatter plots with the parametrizations 
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Figure 4: Evolution of the k parameter as a function of the 
reconstructed energy for experimental data compared with 
pure primary spectra for the angular range 0-24°. The error 
bars assign statistical as well as recostruction uncertainties 
of k. 



dispersion of the k parameter among different bins, which 
include statistical errors, and systematic uncertainty of 
equations 3 in each angular bin. Fig. 4 shows also two thick 
straight lines which are used to classify events into differ- 
ent mass groups. The classification is done for each angu- 
lar bin, independently, to avoid possible systematic effects 
among the bins. The top thick line represents the separation 
between e.p. and intermediate samples and it is defined by 
fitting the k e , p ,(E) = (k S i(E) + k c {E))/2 points which 
are obtained by averaging the values of k for Si and C com- 
ponents. In analogy, the bottom thick line represents the 
separation line between intermediate and e.r. samples and 
it is defined by fitting the k e , r ,{E) = (k c {E) + k He (E))/2 
points which are obtained by averaging the values of k for 
C and He components. In the following, e.p. events will be 
defined as those having a k value higher than the top thick 
line and e.r. events those with k below the bottom thick 
line. The region in between, which is dominated mainly 
by CNO and highly contaminated by Si and He, will be 
defined as the intermediate sample. Such an assignment 
is, therefore, chosen on an event-by-event basis. Naturally, 
the absolute abundances of the events in the three classes 
depend on the location of the straight lines. For that rea- 
son, two thin lines with different slopes bracket the range 
of possible positions of each thick line. They represent the 
uncertainties in defining these energy-dependent selection 
cuts. 

Fig. 5 shows the flux of the three components classified as 
previously explained. The reconstructed spectrum of the 
e.p. sample shows a distinct knee-like feature around 8 x 
10 16 eV. Applying a fit of two power laws to the spectrum 
with a smooth transition [13] results in a statistical signifi- 
cance of 3.5cr that the entire spectrum cannot be fitted with 
a single power law. The change of the spectral slope is A7 
= -0.48 from 7 = -2.76 ± 0.02 to 7 = -3.24 ± 0.05 with the 
break position at logio(E/eV) = 16.92 ± 0.04. Applying 
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Figure 5: Reconstructed energy spectrum of the ex, e.p. 
and intermediate samples as classified in the text together 
with the all-particle spectrum. 



the same function to the all-particle spectrum results in a 
statistical significance of only 2.1a that a fit of two power 
laws is needed to describe the spectrum. Here the change 
of the spectral slope is from 7 = -2.95 ± 0.05 to 7 = -3.24 
± 0.08, but with the break position again at logio(E/eV) = 
16.92 ± 0.10. The spectrum of the e.r. sample is compati- 
ble with a single power law (7 ~ -3.37). 
The error bands in the spectra show that the result is in- 
dependent from the particular choice of the selection-cut 
lines. However, in order to further validate the present 
result, parallel shifts of the cut-line of the e.p. sample 
have been applied. Specifically, by shifting the cut line to 
higher values of k, the e.p. sample is enhanced, and its flux 
diminuishes, while shifting the cut line towards lower val- 
ues of k, the e.p. sample becomes more contaminated by 
e.r. events and the flux increases. Fig. 6 indicates that shift- 
ing up the line cut keeps the knee-like structure unchanged, 
while shifts down tend to smooth out the structure. This 
result confirms that the structure seen in the spectrum is 
caused by the e.p. component, and that the conclusion is 
essentially independent of the particular hadronic interac- 
tion model used in the analysis. 

4 Conclusions 

A first result on the evolution of the k parameter, which 
is sensitive to the evolution of the average mass compo- 
sition has been derived using an event-by-event approach, 
based on QGSJet 11-03 model. It shows that the e.p. com- 
ponent shows a knee-like structure at (8.3 ± 0.8) X 10 16 
eV and the change of slope is about A7 ~ 0.5 through the 
break. The change of slope of A7 ~ 0.3 observed in the 
all-particle spectrum is located in the same energy range 
where the change of slope of the e.p. component occurs, 
indicating that it reflects the knee in the e.p. component. 
The e.r. component has a much steeper spectrum com- 
pared to the e.p. one. In the framework of the QGSJet- 



Figure 6: Effect of the shift of the cut line of the e.p. com- 
ponent on its energy spectrum. 

11-03 model, the e.p. component is the most dominant one 
in the energy range 10 16 - 10 18 eV. This result is in agree- 
ment with those obtained by other analyses of KASCADE- 
Grande data [14, 15, 16]. 
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Abstract: In this work, the shower ratio Y CIC = log N^/ log N cn . between the muon and the charged numbers, both 
corrected by atmospheric attenuation with the Constant Intensity Cut method (CIC), is studied as a parameter to separate 
the KASCADE-Grande data into different mass groups. MC simulations performed with CORSIKA on the framework 
of FLUKA/QGSJet II are employed to obtain the expected Y CIC distributions as a function of the energy for different 
cosmic ray primaries as a basis for the separation. Then the Y CIC parameter is used to divide the KASCADE-Grande 
data into "electron rich" and "electron poor" events. The fraction of events generated by H (Fe) primaries and classified 
as "electron rich" ("electron poor") and its dependence on the primary energy is discussed. Finally, the energy spectra 
of these samples are reconstructed. A knee-like structure is found around 10 17 in the "electron poor" sample associated 
with the heavy component of cosmic rays. 

Keywords: KASCADE-Grande, energy spectrum, cosmic ray composition 



1 Introduction 

In order to understand the origin, nature, propagation and 
acceleration mechanism of galactic cosmic rays, detailed 
measurements with enough statistics on their energy spec- 
trum, composition and arrival directions are needed. The 
task is complicated since at high energies (> 1 PeV), cos- 
mic rays must be studied indirectly by means of exten- 
sive air showers (EAS) detected with Earth-bound exper- 
iments. Information about composition and primary en- 
ergy requires the simultaneous observation of several EAS 



parameters, such as the muon and the electron particle 
contents at ground level. By combining the informations 
from both EAS observables, the KASCADE experiment 
was able to separate for the first time the cosmic ray energy 
spectra of different mass groups in the so called knee en- 
ergy region (around 10 15 eV) [1]. The picture of the galac- 
tic cosmic ray spectrum is still incomplete due to the lack 
of statistics in the interval 10 16 — 10 18 eV, where a knee in 
the heavy component of cosmic rays is predicted by several 
models. To explore this energy regime KASCADE was up- 
graded to KASCADE-Grande, a multidetector-setup with 
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enhanced area also designed to measure with precision the 
charged, muon and electron contents of air showers for cos- 
mic ray studies [2]. In this paper we present an analysis 
based on the ratio between the muon and the charged par- 
ticle numbers in EAS. This ratio is exploited to divide the 
events in two samples: the electron poor and the electron 
rich events. The two samples can be identified with EAS 
generated by primaries belonging to different mass groups 
(i.e. light and heavy elements). The technique is illustrated 
by applying it to the KASCADE-Grande data to reconstruct 
the light and heavy spectra of cosmic rays. 

2 The KASCADE-Grande experiment 

KASCADE-Grande, located in Karlsruhe, Germany (110 
m a.s.l.), is a cosmic ray detector designed to measure EAS 
within the energy interval of 10 16 — 10 18 eV. Two detec- 
tor arrays of different sizes are the main components of the 
experiment [2]. The first one is composed by 37 plastic 
scintillator detectors and covers a surface of 700 x 700 m 2 . 
It is aimed to provide measurements of the charged parti- 
cle number (N c h), core position and arrival direction of air 
showers[2] . The penetrating component (N^) is measured 
with the aid of a smaller array of 200 x 200 m 2 integrated by 
252 e/7 and shielded detectors [3, 2]. The electron shower 
size is obtained from the difference between N c h and Np 
parameters. 

3 Simulations and data selection 

To perform the present analysis, MC simulations must be 
invoked. FLUKA [4] and QGSJet 11-03 [5] were employed 
to describe the hadronic interactions at low and high ener- 
gies, respectively. CORSIKA [6] was used to describe the 
EAS development and GEANT 4 to simulate the response 
of the KASCADE-Grande experiment to the passage of the 
air shower. Both simulated and experimental data are saved 
with identical output formats and analyzed with the same 
reconstruction program. MC data was generated for H, He, 
C, Si and Fe nuclei. Events were sampled from an isotropic 
distribution with spectral index 7 = —2. The simulated 
data was weighted to emulate a 7 = — 3 energy spectrum. 

With the aid of Monte Carlo simulations, quality cuts were 
investigated. They were carefully selected to reduce the 
systematic uncertainties on the muon and charged particle 
numbers. First, the experimental data sample was built out 
of events successfully reconstructed with the KASCADE- 
Grande procedure and collected during stable periods of 
data acquisition . Only EAS detected with more than 35 
active stations were considered in this work. Then bor- 
der effects were avoided by picking events with shower 
cores located inside a central area of 1.37 x 10 5 m 2 in the 
Grande array. Finally, the analysis was restricted to EAS 
with arrival zenith angles lower than 30°, log N c h > 6 
and log Np > 5.1. The value is corrected for sys- 
tematic effects by means of a proper correction function 
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Figure 1 : Estimations of the Y CIC ratio for five different 
primaries shown as a function of the reconstructed energy. 
The results are obtained by means of a full EAS (based on 
the QGSJet 11-03 high energy hadronic interaction model) 
and detector simulation. 

parameterized out of MC simulations. With the above se- 
lection cuts, the effective time of observation is equivalent 
to 1173 days. Full efficiency (> 95%) is achieved for 
log(£/GeV) > 7.4. 

4 The Y CIC ratio 

The KASCADE-Grande (N^, N ch ) data offers a unique 
opportunity to investigate the chemical composition of 
galactic cosmic rays. Several methods can be envisaged 
for this enterprise. The one described in this contribution 
relies on the classification of the EAS on the basis of the 
yCic rat j ^ a parameter defined as 

Y CIC = log N^Oref)/ log N ch (6 ref ). (1) 

In the above formula, the EAS observables, and N c h, 
have been corrected event by event for attenuation effects 
in the atmosphere in such a way that they correspond to 
the shower sizes at a reference zenith angle, 6 re f. The cor- 
rection is achieved by applying the Constant Intensity Cut 
Method (CIC) as described in references [7] and [8]. This 
procedure allows to combine data collected from different 
zenith angles in a model independent way. Here, re f was 
chosen to 22°. 

The Y CIC parameter is a sensitive quantity to the com- 
position of primary cosmic rays. The point can be better 
appreciated by looking at figure 1 . There the mean values 
of Y CIC are plotted as a function of the reconstructed en- 
ergy for five single primaries (see [10] for details about the 
energy estimation), the error bars represent the RMS of the 
yCic distributions. We can notice that at a fixed energy, 
the Y CIC ratio grows with the mass of the primary parti- 
cle. The effect can be understood as a result of the fact that 
the physics of hadronic interactions favors the production 
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Figure 2: Fraction of hydrogen (iron) events classified as electron rich (electron poor) using the cut Y^J. C = 0.84. The 
fractions for protons (stars) and iron nuclei (dots) are plotted against the reconstructed primary energy and the distance of 
the core position to the center of KASCADE. These results are obtained with a full EAS simulation based on the QGSJet 
11-03 interaction model. 



of more secondary charged pions, and hence of muons, for 
heavy nuclei. This implies a lower production of electrons 
and charged particles in general when increasing the mass 
number of the parent nuclei. As seen in figure 1, the Y CIC 
distributions are almost energy independent in the region of 
100% efficiency. Taking advantage of the overall behavior 
of Y CIC and in a first attempt to study the primary com- 
position of cosmic rays in KASCADE-Grande by means of 
this parameter a cut is applied on the data at Y^J. C = 0.84. 
Therefore the events are divided into two different sets for 
a subsequent analysis: the electron rich and electron poor 
groups corresponding to events with Y CIC < 0.84 and 
yCic > 0.84, respectively. From figure 1, it is clear that 
the first set covers the mean Y CIC values for light mass 
elements, here H and He, while the second set contains the 
region for heavy mass nuclei, such as Si and Fe. In this 
way the electron rich and electron poor samples become 
representative of the light and heavy mass groups. 

We remind that the Y CIC distributions were obtained in 
the framework of the QGSJet 11-03 hadronic model. How- 
ever other hadronic interaction models, such as EPOS 1 .99 
[9], were investigated. Through these studies it was ob- 
served that the mean values of Y CIC still do not depend on 
the primary energy, while their absolute values are modi- 
fied. Thus, in the case of the EPOS 1 .99 hadronic interac- 
tion model, the optimal separation between light and heavy 
mass groups is Y^ c = 0.86. 

The separation into different mass groups employing the 
yCic cut j s not com pi e tely clean due to the size of the 
fluctuations inside each energy bin. Using MC simulations 
we have evaluated, for each element, the fraction of events 
selected as electron rich or electron poor. For simplicity, 
only the results obtained for hydrogen and iron nuclei are 
shown in figure 2 as a function of the reconstructed en- 
ergy (left panel) and the core position (right panel). For 
energies greater than \og(E/GeV) > 7.4 (i.e. 100% ef- 



ficiency) the fraction of events properly classified through 
the Y CIC ratio is almost energy independent and greater 
than 80%. This fraction is somewhat higher for iron pri- 
maries because of the lower shower fluctuations for heavy 
nuclei. It is not surprising to find out that their correspond- 
ing curves are approximately constant (in the energy range 
of full detection efficiency) over the whole effective area 
since the selection cuts on the EAS parameters were explic- 
itly chosen to guarantee that measured shower parameters 
do not depend on the experimental conditions, such as the 
core distance from the KASCADE array center (i.e. mean 
distance of the muon detectors from the shower core). 

These results show that possible distortions introduced in 
the reconstructed energy spectra of the different event sam- 
ples, separated using the Y^ c value, are minimized. In 
the next section, the method is used to get insight into the 
composition of the primary cosmic ray flux. 

5 Applications: Energy spectra for two mass 
groups 

Once the experimental data is divided into electron poor 
and electron rich groups using the above procedure above, 
the primary energy of each event is estimated following 
[10]. After the energy derivation, the energy spectrum is 
reconstructed for each electron group leading to the energy 
spectra for the light and heavy components of the cosmic 
ray primary flux, which are presented in figure 3. The all- 
particle spectrum is also shown, and it is estimated as the 
direct sum of the two mass group spectra here estimated. 
The plots are not corrected for migration effects. 

From figure 3 an outstanding feature is immediately ap- 
preciated in the heavy component of cosmic rays which 
corresponds to a knee-like structure. This feature also ap- 
pears in more detailed analyses of the KASCADE-Grande 
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Figure 3: Light and heavy energy spectra derived from the 
KASCADE-Grande (N^Ncf,) data using the Y C ic ratio 
as a mass estimator. The all-particle cosmic ray energy 
spectrum (squares) derived by adding the individual mass 
group fluxes is also shown. 

data [12, 13, 14]. A fit with a broken power-law spec- 
trum to the electron poor spectrum indicates a change of 
slope A 7 = -0.47 at log 10 (E/GeV) = 7.77 ± 0.06 from 
7 = -2.72 ± 0.02 to 7 = -3.19 ± 0.02. As seen in figure 
4, the knee-like feature is a structure which does not de- 
pend on the position of the Y~Jj. , i.e., it is inherent to the 
cosmic ray data. 

In figure 3, it can be noticed that the all-particle energy 
spectrum shows also a break. From a fit with a broken 
power-law spectrum, the position of the break is found 
around log 10 (S/GeV) = 7.92 ± 0.10, but the change of 
slope is smaller for this case: A7 = —0.29, from 7 = 
-2.95 ± 0.05 to 7 = -3.24 ± 0.08. It is seen that this fea- 
ture arises as a consequence of the knee in the heavy com- 
ponent, but is less pronounced due to the contribution of 
the light group, which is not zero for log(i?/GeV) > 7.4. 

6 Conclusions 

In this contribution, the potential of the Y CIC ratio as a 
mass sensitive parameter was studied. It was shown that 
the classification of EAS data into electron rich and elec- 
tron poor events by applying the cut Y~^. c = 0.84 can lead 
to the separation of the light and heavy mass groups in the 
cosmic ray flux. The fraction of events generated by hydro- 
gen (iron) primaries and classified as electron rich (elec- 
tron poor) does not depend on the primary energy and (in 
the frame of the QGSJet 11-03 interaction model) is greater 
than 80%. Properly choosing the selection cuts on the EAS 
parameters the separation efficiency can be expected to be 
constant with energy, core position and zenith angle avoid- 
ing the introduction of distortions to the corresponding en- 
ergy spectra. By applying the Y CIC method, the energy 
spectra of the light and heavy mass groups were obtained. 



Figure 4: Energy spectra for the electron poor sample ob- 
tained for different Y CIC cuts. 

The results show that the heavy component presents a knee 
around log 10 (£'/GeV) = 7.77 ± 0.06 and that the con- 
tribution to the all-particle energy spectrum from the light 
component is different from zero. 
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Abstract: The KASCADE-Grande experiment, located at Karlsruhe Institute of Technology, is a multi-component Ex- 
tensive Air Shower (EAS) detector studying primary cosmic rays in the 10 16 — 10 18 eV energy range. In this contribution 
a measurement of the cosmic ray chemical composition, based on the comparison of the experimental distributions of 
the ratio between the EAS muon size and the total charged size (N M /N c h) with those expected from a complete EAS 
simulation based on the QGSJet 11-03 interaction model, is presented. It has already been shown that, in the frame of 
this interaction model, the detector performances allow to separate three different mass groups: light, intermediate and 
heavy. With the employed technique, the relative abundances of the three mass groups are derived and their evolution as 
a function of the charged particle size is observed. The differential spectra of the three mass groups are then inferred and 
discussed 

Keywords: Composition, KASCADE-Grande 



1 Introduction 

The energy range from 10 16 to 10 18 eV is very important 
to deeply investigate the details of the knee of the primary 
energy spectrum. In the last ten years different experiments 
(e.g KASCADE [1] and EAS-TOP [2]) have shown that 
this feature can be attributed to the light component of the 
primaries (i.e. H or He). A definitive confirmation of the 
astrophysical origin of the knee can only be obtained by 



the detection of the foreseen change of slope of the heavy 
component at a primary energy ~ 10 17 eV. 

In this contribution we present a study of the primary cos- 
mic ray chemical composition in the 10 16 — 10 18 eV en- 
ergy range performed comparing the KASCADE-Grande 
experiment [3] data with the results of a full Extensive Air 
Shower simulation based on the CORSIKA [4] code and 
on the QGSJet 11-03 [5] interaction model. 

The KASCADE-Grande experiment is located at the Cam- 
pus North of the Karlsruhe Institute of Technolgy, 110 
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m a.s.l. It consists of an array of 37 plastic scintilla- 
tor modules 10 m 2 each (Grande) spread over an area of 
700x700 m 2 , working jointly with the co-located and for- 
merly present KASCADE experiment [6], consisting of 
252 electron and muon scintillation detectors placed over a 
200x200 m 2 area. All events triggering a sevenfold coinci- 
dence of Grande array detectors, arranged on an hexagonal 
grid (mean side lenght of 130 m), are reconstructed. The 
charged particle size N c /j is determined by the Grande array 
and the muon number N M by the KASCADE array. Both 
(N c h) and (N M ) are measured with an accuracy < 15%; a 
detailed description of the reconstruction procedure can be 
found in [3]. 

2 The Chi Square Method 

The present analysis is performed, to minimize effects of 
a possible incorrect description (in the simulation) of the 
EAS evolution in atmosphere, on vertical events (0° < 
9 < 24°). Data are divided in twelve N c /j intervals (from 
Log N c h=6.0 to Log N c /,=8.0) and the N^/N c /j experimen- 
tal distributions are fitted with a linear combination of those 
obtained by a complete simulation of events generated with 
a power law distribution in the 10 15 — 10 18 eV energy range 
for five different primaries (H, He, C, Si, Fe). We have al- 
ready shown that the KASCADE-Grande experiment per- 
fomances allow us to separate, with such algorithm, events 
into three samples originated by different mass groups [7]. 
The combination giving the best fits (examples for two of 
the twelve N c h intervals are shown in figures 1 and 2) is 
the one composed by the following mass groups: light, 
100% hydrogen; intermediate, 50% helium and 50% car- 
bon; heavy, 50% silicon and 50% iron. The linear com- 
bination of the simulated mass groups is expressed as fol- 
lows: 

(1) 



-^sjm(^) ^ ^ ^3 fsim, j 



where F S i m {i) is the total theoretical fraction of simulated 
events falling in the histogram channel i, f s i m ,j(i) is the 
fraction for the single component j (j = 1,2,3), is the 
sum over the different components and atj is the fit param- 
eter representing the relative abundance of the component 
j. The fit parameters fulfill the conditions 



and 



0. < oij < 1. 



J2 a i = L 



(2) 



(3) 



The fit is performed minimizing the following Chi Square 
function: 



E 



(4) 



where F exp (i) is the fraction of experimental events falling 
in the histogram channel i and a(i) is the error on the the- 
oretical 
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Figure 1: The experimental N^/Nch distribution (solid 
line) measured in the 6.24 < Log(N c h) < 6.36 interval 
for vertical events (0° < 8 < 24°) fitted by a combina- 
tion of Hydrogen (thin dashed line), Helium + Carbon in 
a 50% mixture (dot dashed line) and Silicon + Iron in a 
50% mixture (thick dashed line). The triangles show the 
fit result. The experimental plot is normalized to 1 and ev- 
ery simulated component is normalized to its relative abun- 
dance. The tails of the experimental distribution are treated 
summing the events (the corresponding value is shown by 
a star) in a single bin, in order to have at least 5 events 
in each of them. In the shown example, the results for 
the abundances and the fit are: aught = 0.14 ± 0.02; 



Q-intermediate = 0.59 ± 0.02; OLheavy = 0.27 ± 0.01. 

The chi square and the cumulative function values are: 
xl/y = 56.52/46 = 1.23, P( X 2 > xl) = 0-14 



expression (1). A fit is validated and the relative abun- 
dances of the mass groups are accepted if the cumulative 
function, P(\ 2 > xl), is included in the 0.05 < P < 0.95 
interval. 



The choice of sampling events in bins of the charged par- 
ticle size is made in order to have a selection depending 
only on the performances of the experiment and not being 
influenced by an EAS simulation (as could be in the case 
of reconstructed energy bins). 

The behavior of the relative abundances versus the charged 
particle size is shown in figure 3. The light and intermedi- 
ate components show a trend with bigger fluctuations with 
respect to the heavy mass group. The relative abundance of 
the light mass group is almost independent on the charged 
particle shower size, while the heavy one shows a sudden 
change at \ogN c h <~ 6.8. The intermediate mass group 
seems to be more abundant near (in N c h bins) the thresh- 
old, decreasing for higher N c h and becoming nearly similar 
to those of light elements. 

3 Spectra of the single mass groups 

From the number of events measured in each A \og{N c h) 
interval (N exp (k),k = 1, 12) and having determined the 
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Figure 2: The experimental N^/N c h distribution (solid 
line) measured in the 7.15 < Log(N c h) < 7.24 interval 
for vertical events (0° < 6 < 24°) fitted by a combination 
of light (thin dashed line), intermediate (dot dashed line) 
and heavy (thick dashed line). The triangles show the fit re- 
sults. In the shown example, the results for the abundances 
and the fit are: au gh t = 0.16 ± 0.02; a int er m ediate = 
0.30 ± 0.04; a heavy = 0.52 ± 0.03. The chi square and the 
cumulative function values are: x\l v = 9.66/14 = 0.69, 
P(X 2 > Xl) = 0.79 
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Figure 3: Relative abundances of the three mass groups 
versus the charged particle size N c h measured for vertical 
events. 





light 


intermediate 


heavy 


a 


0.99 ± 0.01 


0.93 ± 0.004 


0.91 ±0.01 


b 


0.75 ± 0.06 


1.30 ±0.02 


1.61 ±0.03 



Table I: a, b parameters of the power law correlating, for 
each mass group, the charged particle size N c /j and the pri- 
mary energy E 



relations we convert the N c /j intervals into the correspond- 
ing energy bins (that are thus not the same for the three 
mass groups) and from the measured fluxes we derive the 
differential flux at the central energy of each interval. 

To verify if, through the described algorithm, we are able 
to reproduce the spectra of the single mass groups we have 
performed the same analysis on test spectra: half of the 
simulated data set is used as fake experimental data while 
the other half is used as reference for the fitting procedure. 
With the goal of identifying possible spectral distortions 
artificially introduced by the analysis we have chosen the 
test spectra with a slope 7 — — 3 and without breaks for 
all elements. The results are shown in figure 4: the heavy 
mass group shows good agreement between test and recon- 
structed spectra in the whole energy range. The light and 
intermediate mass groups are in good agreement for en- 
ergies greater than 4 x 10 16 eV while near the threshold 
the reconstructed spectra are lower than the test ones. This 
result is not unexpected as shower and experimental fluc- 
tuations (that are bigger for low masses and originate an 
event migration from bin to bin) are not yet corrected for 
in the conversion from N c h to energy. This procedure is 
under development, its results will be shown soon. It is 
worthwhile to point out that the uncorrected reconstructed 
test spectra show, for all mass groups, no artificial breaks. 

Another hint of the correct behavior of the described analy- 
sis algorithm can be derived from figure 5 showing the rela- 
tive abundances reconstructed for the previously described 
test spectra. These abundances show no dependence on 
N c h in the whole range, even if fluctuations (mainly of the 
light and intermediate mass groups) around the expected 
value can be seen. 



relative abundance of each mass group (J) we can calcu- 
late the number of events that are originated by primaries 
belonging to the j th mass group (Nj(k) = ay • N exp (k)). 
The three mass group fluxes in each charged particle size 
bin are then derived: 

<f>j(k) = iV i (fc)/(ArJA<)(m" 2 sr" 1 s" 1 ) (5) 

where A = 1.52 x 10 4 m 2 , = 0.52sr and At = 1.1 x 
10 8 s. 

From a full EAS (based on the QGSJet 11-03 interaction 
model) and detector simulation we derive, for each mass 
group, a power law relation (Eq = 10 b • N^ h ) between the 
charged particle size and the primary energy (the values of 
the a and b parameters are reported in table 1). With these 



4 Conclusions 

The performances of the algorithm aiming at the measure- 
ment of the evolution of the primary chemical composition 
of cosmic rays are described. We have shown how, with 
the KASCADE-Grande experiment performances, we are 
able to separate three mass groups and that all of them are 
needed to reconstruct the measured N^/N c h distibutions. 
The measured chemical composition gets heavier as the 
shower size N c h increases as can be inferred from the evo- 
lution of the relative abundances of the three mass groups 
with N ch . 

The measured spectra of the three mass groups are shown 
in figure 6 as obtained in the frame of the QGSJet 11-03 
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Figure 4: Comparison between test and reconstructed spec- 
tra for the three mass groups. In the fake spectrum all mass 
groups have equal abundance, in the plot the fluxes of light 
and heavy mass groups are multiplied by arbitrary factors 
for clarity. 
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Figure 6: Differential energy spectra of the three mass 
groups: heavy (Si+Fe, dots), intermediate (He+C, trian- 
gles) and light (H, squares) elements. 

All shown spectra do not yet take into account effects of 
event-to event shower fluctations at the energy assignment: 
a procedure to unfold them is under development. 

It is important to emphasize that the present results, and 
in particular the relative abundances of the mass groups, 
heavily depend on the choice of QGS Jet 11-03 as interaction 
model. The entire analysis will be, in near future, repeated 
on basis of a complete EAS simulation based on a different 
interaction models, i.e. EPOS [10]. 
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Figure 5: Reconstructed relative abundances of the three 
mass groups versus the charged particle size N c h obtained 
analyzing test spectra with constant and equal slope for all 
elements. 

interaction model. The spectrum of the heavy mass group 
shows a distinct change of slope and cannot be described by 
a single power law. This reconstructed spectrum is fitted 
with a function [8] describing a spectral shape with two 
different slopes interconnected by a smooth knee at energy 
Ek- The results, shown by the solid line in figure 6, are: 
7i = 2.67 ± 0.02, 72 = 3.29 ± 0.02 and \ogE k {eV) = 
17.79 ± 0.04. The statistical significance of the change of 
slope is ~ 4cr. At a similar energy a (weaker) steepening 
is also observed in the all particle spectrum measured by 
KASCADE-Grande [9]. 

In contrast both the intermediate and light mass groups 
spectra are descriable by a single power law due to the low 
event numbers. These spectra are dominated by statistical 
fluctuations, and no significant conclusion can thus be de- 
rived at the moment. 
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Abstract: The KASCADE-Grande experiment, located at KIT-Karlsruhe, Germany, consists of a large scintillator array 
for measurements of charged particles, N c h , and of an array of shielded scintillation counters used for muon counting, iV M . 
KASCADE-Grande is optimized for cosmic ray measurements in the primary energy range 10 16 eV to 10 18 eV, thereby 
enabling the verification of a possible second knee expected at approximately 10 17 eV. Exploring the composition in this 
energy range is of fundamental importance for understanding the transition from galactic to extragalactic cosmic rays. 
Following earlier studies of elemental spectra reconstructed in the first knee energy range from KASCADE data, we shall 
now extend these measurements to beyond 10 17 eV. By analyzing the two-dimensional shower size spectrum N cn , vs. 
Np, we reconstruct the energy spectra of different mass groups by means of unfolding methods. The procedure and its 
results, which yield a strong indication for a kink in the iron spectrum at around 80 PeV, will be presented. 

Keywords: Cosmic ray, energy spectrum, composition, knee, iron knee, KASCADE-Grande 



1 Introduction 

The spectrum of cosmic rays follows a power law over 
many orders of magnitude in energy, overall appearing 
rather featureless. However, there are a few structures ob- 
servable. In 1958 Kulikov and Khristiansen [1] discovered 
a distinct steepening in the spectrum at around 10 15 eV. 
Three years later Peters [2] concluded that the position of 
this kink, also called the cosmic ray "knee", will depend 
on the atomic number of the cosmic ray particles if their 
acceleration is correlated to magnetic fields. Round about 



half a century later, the KASCADE experiment [3] clarified 
that this change in spectral index is caused by a decrease of 
the so far dominating light 1 component of cosmic rays [4]. 
This result was achieved by means of an unfolding anal- 
ysis disentangling the manifold convoluted energy spectra 
of five mass groups from the measured two-dimensional 
shower size distribution of electrons and muons at obser- 
vation level. Based on the high energy interaction model 
QGSJET 01 [5] it was shown, that the kink in the all- 

1. The description "light" refers to the atomic mass of the cos- 
mic ray particles, which are primarily nuclei. 
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Figure 1 : Measured shower size distribution. 

particle spectrum at around 5 x 10 15 eV corresponds to 
a knee observed in the hydrogen flux. 
Nowadays, there are numerous theories about the origin 
and acceleration of cosmic rays. Concerning the knee po- 
sition, some of them predict in contrast to the rigidity de- 
pendence considered by Peters a correlation with the mass 
of the particles. Hence, it is of great interest to verify 
whether also the spectra of heavy components exhibit anal- 
ogous structures and if, at what energies. The KASCADE- 
Grande experiment [6] extends the accessible energy range 
of KASCADE to higher energies up to 10 18 eV and allows 
by this to investigate the cosmic rays composition at re- 
gions where the possible, so-called iron knee is expected. 
The determination of this "second" knee enables the val- 
idation of the various theoretical models. Following this 
purpose, the KASCADE-Grande measurements have been 
analyzed similar to the afore-mentioned studies [4] of the 
KASCADE data. The applied unfolding method will be 
outlined in the next section. Thereafter, the uncertainties of 
the analysis and the resulting elemental energy spectra will 
be shown and studied. A more comprehensive description 
can be found in [7]. 

2 Outline of the analysis 

The analysis' objective is to compute the energy spectra 
of five 2 cosmic ray mass groups, represented by hydro- 
gen (H), helium (He), carbon (C), silicon (Si) and iron 
(Fe) nuclei, from 10 16 eV beyond primary energies of 
10 17 eV. The two-dimensional shower size spectrum \gN c h 
vs. lgiV^ of charged particles and muons measured with 
KASCADE-Grande is used as starting point for the un- 
folding analysis (Fig. 1). All measured 3 air showers have 
to pass certain quality cuts to ensure good reconstructed 
shower sizes. Particularly, only air showers with zenith an- 
gles less than 18° are used exhibiting at least 10 6 charged 
particles and 10 5 muons. The measurement time covers ap- 
proximately 1 318 days resulting in 78 000 accepted events 
passing all quality cuts. This corresponds to an exposure of 
164 709 m 2 sr yr. 

The convolution of the sought-after differential fluxes 
dJ n /d\gE of the primary cosmic ray nuclei n into the 
measured number of showers Ni contributing to the cell i 
of shower size plane, and thus to the content of this specific 



charged particle and muon number bin (\g(N c h), lg(iV A1 )) i , 
can be described by an integral equation: 

N i = 12 J J J J jte Pn dlgE cos9dA dn dt > 

™ =1 T m Otat At E 

(1) 

with 

Pn=Pn((lgN ch ,lgN tl ) i \ lgE) . 

One has to sum over all N n elements contributing to the 
all-particle cosmic ray spectrum, in this analysis the five 
representative primaries. T m is the measurement time, Qtot 
the total solid angle accessible for the experiment and used 
for the analysis, and At the chosen fiducial area. The term 
p n represents the conditional probability to reconstruct a 
certain combination of charged particle and muon number, 
respectively to get an entry in the cell (lg(AT c ^), ]g(JV^)) ., 
if the air shower inducing particle was of the type n and 
had an energy of E. More precisely, p n itself is a convolu- 
tion combining the intrinsic shower fluctuations occurring 
whilst the air shower development, the detection and recon- 
struction efficiency as well as the properties of the observ- 
ables' reconstruction process. The cosine term in cos 9 dA 
accomplishes the transformation from the horizontal sur- 
face element to the effective detection area. 
Equation (1) can mathematically be understood as a sys- 
tem of coupled integral equations referred to as Fredholm 
integral equation of first kind. There are various methods 
to solve such an integral equation, albeit a resolvability of- 
ten doesn't per se imply uniqueness. In some preliminary 
tests it was found, that the unfolding algorithm of Gold 
[8] yields appropriate and robust solutions. It is a iterative 
procedure and de facto related to a minimization of a chi- 
square function. For countercheck purposes all results are 
validated by means of two additional algorithms, an also 
iterative method applying Bayes' theorem [9] performing 
very stable, too, and an regularized unfolding based on a 
combination of the least-squares method with the principle 
of reduced cross-entropy [10], that yields slightly poorer 
results. All these solution strategies have in common that 
the response 4 function p n of Eq.(l) has to be known a 
priori. It is parametrized based on Monte Carlo simula- 
tions. The air shower development is simulated by means 
of CORSIKA [11] 6.307 based on the interaction models 
QGSJET-II-02 [12] and FLUKA 2002.4 [13]. The experi- 
ment's response is simulated using CRES 5 1.16/07, which 
bases on GEANT 3.21 [14] detector description and simu- 
lation tool. 

2. Due to effects of limited resolution not any number of mass 
groups can be treated. 

3. Most of the cuts, e.g. the chosen zenith angle range, are also 
applied to the simulated air showers. 

4. Also named kernel or transfer function; and more precisely 
it is rather a matrix than a simple function. 

5. Cosmic Ray Event Simulation, a program package devel- 
oped for the KASCADE [3] detector simulation. 
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3 Error analysis 

The determination of the elemental energy spectra will be 
subjected to influences of different error sources. They can 
roughly be classified in two categories: uncertainties in- 
duced, or at least appearing whilst the deconvolution pro- 
cess and those embedded in the computed response func- 
tion caused by the limited Monte Carlo statistics. 

3.1 Uncertainties whilst the deconvolution 

Firstly, the used data set is only a small sample based on 
a limited exposure, and hence suffering from statistical un- 
certainties. They are propagated through the unfolding al- 
gorithm and affect the quality of the solution. Furthermore, 
the used deconvolution method itself can introduce a sys- 
tematic bias. The influences of both sources can be evalu- 
ated by means of a frequentist approach. Assuming appro- 
priate 6 spectral indices some trial elemental energy spectra 
are specified based on which a test data sample can be gen- 
erated using Eq.(l). Subsequently, these data samples are 
unfolded. Since the true solution is a priori known, the de- 
convolution result can be compared to it to reveal statistical 
fluctuations induced by the limited measurement time and a 
possible systematic bias induced by the unfolding method. 

3.2 Influences of limited Monte Carlo statistics 

The amount of simulated air showers is strongly limited 
by reason of computing time. Due to the limited Monte 
Carlo statistics, the computation of the response function, 
i.e. the parametrization of the intrinsic shower fluctuations 
as well as of the detector properties, will only be possi- 
ble under certain uncertainties resulting in an systematic 
error of the finally unfolded solution. In Fig. 2 exemplar- 
ily the simulated charged particle number distribution in 
case of hydrogen induced air showers with primary energy 
of 2 x 10 15 eV is shown. A scattering around the used 
parametrization ("normal") can be observed. This statisti- 
cal uncertainty will be treated conservatively: Considering 
the computed fit parameters and their errors some new sets 
of parameters are calculated by means of a random genera- 
tor. Based on each set, new response functions can be com- 
puted and used to unfold the data. Comparing the results 
reveals the caused systematic uncertainty in the solution. 
The distributions' tails have to be inspected in more detail. 
Because of the very low statistics, the tails can vary within a 
certain range without worsen the fit result. In particular the 
right tail describing the fluctuations in direction to higher 
energies can have an important impact on the unfolded so- 
lution due to the steeply falling flux of cosmic rays. The 
systematic influence of the tails will conservatively be esti- 
mated by computing two additional response functions as- 
suming in contrast to the standard case, in agreement with 
the statistical uncertainties, either a very fast decreasing or 
an elongated tail (cf. Fig. 2). Using both for a deconvo- 
lution and comparing the results yields the maximal sys- 
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Figure 2: The simulated charged particle number distribu- 
tion of 2 x 10 15 eV hydrogen induced air showers with 
angles of incidence less than 18°. The distribution is fitted 
based on different approaches (see text). 

tematic error range caused by the uncertainty in the tails 
description. 



4 Results and conclusion 

In Fig. 3 left panel the unfolded differential energy spec- 
tra of hydrogen and iron as well as the sum 7 of the three 
single spectra of helium, carbon, and silicon are shown, 
representing respectively the fluxes of the light, heavy, and 
intermediate mass groups of cosmic rays. In addition, all 
five unfolded spectra are summed up to the all-particle flux, 
which is compatible to the results of other experiments 
(right panel) and agrees very well with the KASCADE- 
Grande all-particle spectrum published in [15]. The shaded 
band indicates the methodical uncertainties while the error 
bars represent the statistical error originating from the lim- 
ited measurement time. With increasing energy the heavy 
component gets the dominant contributor to the cosmic ray 
composition. This agrees with the results of KASCADE 
[4] where a reduction of the light component beyond the 
first knee was found. 

Both in the all-particle and the iron spectrum there is by 
eye a slight bending discernible at around 10 17 eV. How- 
ever, the change in the all-particle spectrum reveals not to 
be significant. In this context, one should keep in mind 
that this spectrum is the sum of all five elemental spectra 
unfolded separately, and by this is affected by their uncer- 
tainties. For the determination of the all-particle spectrum 
with KASCADE-Grande, there are more precise methods 
available, e.g. that one introduced in [15] stating a high 
significance for a change in the spectral index of the all- 

6. Spectral indices close to those estimated by KASCADE are 
used in order to have realistic ones. However, some unlikely spec- 
tra are tested, too. 

7. The intermediate component was combined because of the 
poorer results suffering from low statistics and making the plot 
confusing without giving further insights. 
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Figure 3: Depicted are the unfolded energy spectra for H and Fe, a combined spectrum for He, C and Si, as well as the 
all-particle spectrum (left panel). The all-particle spectrum conforms well to those of other experiments (right panel). 



particle spectrum at around 1 x 10 17 eV. 
In order to judge the possible structures in the unfolded 
iron spectrum, it was fitted preliminarily by a single power 
law. However, the resulting chi-square probability for such 
a featureless single power law was below 1% (x 2 /ndf = 
18.9/7). In Fig. 4, the residual flux between the iron spec- 
trum shown in Fig. 3, left panel, and such a spectrum that 
was derived by a single power law fit is depicted in order 
to emphasize the deviations between the single power law 
and the unfolded spectrum. Additionally, the iron spectrum 
is now fitted by a double power law: 



dJ(E) 
dlgE 



po x E P2 x 1 



P4\ (P3-P2)/P4 



(2) 



where p\ = \g(E-k ncc /GeW) = 7.9 ± 0.1 corresponds to 
the knee position, while p2 = —2.62 ± 0.02 and P3 = 
—3.7 ± 0.4 are the spectral indices below and above the 
knee. The sharpness of the knee structure is encoded in 
P4 = 7.0 and was fixed without worsening the fit's qual- 
ity, while po is a free normalization parameter. This fit de- 
scribes the spectrum significantly better (chi-square proba- 
bility at around 30% with x 2 /ndf — 6.2/5), giving strong 
indications for a kink in the iron flux at around 80 PeV. 
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Comparing the position of this potential iron knee to that 
for hydrogen (at around 2 PeV to 4 PeV) gives indications 
for a scaling of the knee positions with the charge of the 
nuclei rather than with their atomic mass number. This 
would encourage the cosmic ray acceleration models based 
on magnetic fields. 

To summarize, there is a strong indication for a kink in the 
iron-like spectrum at around 80 PeV as well as for a depen- 
dence of the cosmic ray acceleration process on the charge 
of the nuclei, both on the premise that especially the used 
model QGS JET-II-02 describes the physics of hadronic in- 
teractions at these energies with a high level of reliability. 
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Figure 4: Residual flux between the iron spectrum shown 
in Fig. 3 and a spectrum that was derived by a single power 
law fit to that iron spectrum. Additionally, the iron spec- 
trum is now fitted by a double power law. 



8. And assuming that the mass groups represented by H and Fe 
actually consist only, or at least primarily, of that two primaries. 
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Abstract: We present a method to reconstruct the primary energy spectrum of cosmic rays from the charged particle 
densities recorded with the KASCADE-Grande detector. The KASCADE-Grande is hosted by the Karlsruhe Institute 
for Technology (KIT), Germany and is operated by an international collaboration. It has been shown that the charged 
particle density becomes independent of the primary mass at certain fixed distances from the shower axis and that it 
can be used as an estimator of the primary energy. Such distance is a characteristic of the detector and in the case of 
the KASCADE-Grande experiment it was shown that it is 500 m, hence the notation S(500). A relation is established 
by means of simulations between the primary energy of cosmic rays and the S(500). We account for the attenuation of 
inclined showers by applying the CIC method. By using the simulation derived calibration we build the primary energy 
spectrum from the recorded S(500). Several sources of systematic uncertainties are identified and their contribution to 
the final result is evaluated. Additionally we apply an unfolding to account for statistical fluctuations. The features of the 
obtained spectrum are discussed in relation to the result of another independent technique that is applied at KASCADE- 
Grande. 

Keywords: KASCADE-Grande, primary energy spectrum, S(500) 



1 Introduction 

Previous investigations have shown that the charged parti- 
cle density in air showers becomes independent of the pri- 
mary mass at large but fixed distances from the shower axis 
and that it can be used as an estimator for the primary en- 
ergy [!]■ Such a distance is characteristic for a given ex- 



periment. Based on this property a method was derived to 
reconstruct the primary energy spectrum from the particu- 
lar value of the charged particle density, observed at such 
specific radial distances. While in the AGASA experiment 
the technique was applied for a distance of 600 m to the 
shower axis [2], in the case of the KASCADE-Grande ar- 
ray, detailed simulations [4] have shown that the particular 
distance for which this effect takes place is about 500 m 
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Figure 1: Simulations show that, for the case of the 
KASCADE-Grande experimental layout, the particle den- 
sity becomes independent of the primary mass at around 
500 m distance from the shower core; this plot shows av- 
eraged simulated lateral distributions for different primary 
types with equal energy. 



(see Fig. 1 and Fig. 4). Hence the notation S(500) for the 
charged particle density at 500 m distance from the shower 
core. The distance is measured in a plane normal to the 
shower axis and containing the shower core. The data 
recorded in the detector plane is projected on the normal 
plane taking into account the attenuation effects character- 
istic to inclined events. 

The study has been performed for both simulated (Fig. 5) 
and experimental (Fig. 6) events, using identical recon- 
struction procedures [5]. The CORSIKA Monte Carlo EAS 
simulation tool [6] is used to simulate air showers, with the 
QGSJETII model embedded for high energy interactions 
[7]. 

The energy deposits of particles in the KASCADE-Grande 
detector stations are recorded along with the associated 
temporal information (arrival times of particles). Using 
appropriate Lateral Energy Correction Functions (LECF), 
the energy deposits are converted into particle densities. 
The LECF functions are dependent on the shower zenith 
angle [8] and on the position of the station around the 
shower core (i.e. the LECF are dependent on the angle 
of incidence of particles in detectors). For every recorded 
event, a Linsley [9] Lateral Density Function (LDF) is used 
to evaluate the particle density at the radial range of inter- 
est, 500 m. 

The described reconstruction is performed independently 
from the standard reconstruction applied at KASCADE- 
Grande, based on the N C ^-N M approach and described in 
[10]. 
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Figure 2: S(500) reconstruction efficiency for different 
zenith angular ranges and for the entire shower sample (all 
quality cuts applied); the reconstruction efficiency exceeds 
95% atlogio[£yGeV]>7.5 

2 KASCADE-Grande 

In its general structure and operation, the KASCADE- 
Grande array [3] resembles other cosmic ray experiments. 
Hosted by the Karlsruhe Institute of Technology (KIT), 
Campus North, Germany at 110 m a.s.l and operated by 
an international collaboration, it is composed of many 
detector stations that are distributed over a wide area and 
it has been designed to record charged particle densities, 
without disentangling the type of secondary that is inter- 
acting in the sensitive medium. The shape of the array is 
rectangular with a length of ss700 m. The lateral particle 
density distribution is subsequently inferred by adjusting 
the data registered with the detector stations to an a-priori 
assumed lateral distribution function. 
Historically, the KASCADE-Grande detector array is 
an extension of a smaller array (the KASCADE array, 
operated since 1996). KASCADE was designed to record 
air showers initiated by primaries with energies in the 
10 14 -10 17 eV range (including the knee range). The 
extension of the original KASCADE array was guided by 
the intention to extend the energy range for efficient EAS 
detection to 10 16 -10 18 eV (Fig. 2). This extended energy 
range provides various interesting aspects: the expected 
transition from galactic to extragalactic cosmic rays and, 
in particular the question whether there exists a "second 
knee" in the energy spectrum. 



3 The constant intensity cut method (CIC) 

For a given event sample, an EAS observable could have 
different values for events induced by identical primaries 
(Eq,Ao), but arriving from different zenith angles (due to 
EAS attenuation through the atmosphere). This is also the 
case for the S(500). To reach the detectors at ground level, 
an inclined event will propagate along an extended path 
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Figure 3: Integral S(500) spectra; the horizontal line is a 
constant intensity cut at an arbitrarily chosen intensity; by 
assuming an exponential attenuation pattern the attenuation 
length of S(500) was evaluated to A(500)=754 ± 8 g-cm- 2 . 

through the atmosphere compared to a vertical shower. The 
recorded particle densities at a given radial range for the 
two events will be different. One has to correct for this 
effect before performing an analysis simultaneously on all 
EAS events. This is achieved by applying the Constant 
Intensity Cut (CIC) method (Fig. 3) [11]. The method is 
based on the assumption that for a given minimum primary 
energy threshold we should record the same flux of events 
(i.e. primaries) from all zenith angles. Therefore we di- 
vide our shower sample in sub-samples corresponding to 
different zenith angular bins that subtend equal solid an- 
gles. Then we look in the integral S(500) spectra of each 
sub-sample at a given arbitrary intensity. We establish a 
correlation between the S(500) corresponding to the given 
intensity in each spectrum and the corresponding zenith an- 
gle. This correlation we use as a correction function. All 
reconstructed S(500) values are corrected for attenuation 
by bringing them to the value they would have at a chosen 
reference angle. For the present study the reference angle is 
considered to be 21°, since the zenith angular distribution 
for the recorded EAS sample peaks at this value. The CIC 
correction is derived entirely from recorded experimental 
data and is independent from simulated studies. 

4 Conversion to energy 

For the experimental EAS sample, the total time of acqui- 
sition was 1 173 days for a 500 x 600 m 2 fiducial area. The 
same quality cuts were used for both simulated and experi- 
mental events. Only those events are accepted for which 
the zenith angle is below 30°, the reconstructed shower 
core is positioned inside the detector array and not too 
close to the border, and the event is triggered by more than 
24 Grande stations (to ensure that density data is recorded 
close to the 500 m radial range, i.e. the shower is large 
enough). A good quality of the fit to the Linsley distribu- 
tion is a further important criterion. 

A calibration of the primary energy Eq with S(500) was 
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Figure 4: The dependence of S(500) on the primary energy 
Eo for two different primaries (showers in fairly equal pro- 
portions for the two masses); the box-errors are the errors 
on the spread; the errors on the mean are represented with 
bars and are dot-sized; straight lines represent power law 
fits. 

derived from simulations (see Fig. 4). For the systematic 
contribution to the total error, several sources of system- 
atic uncertainties have been identified: the spectral index 
of the simulated shower sample (which is different from 
the true one) is acting as a source of systematic uncertainty 
(<1% contribution), the S(500)-E calibration (<1% con- 
tribution), the CIC method (<1% contribution), the statis- 
tical fluctuations in the simulated shower sample (7%) and 
the choise of a certain reference angle at which to perform 
the S(500) attenuation correction (7% contribution). 
The energy resolution has also been evaluated from sim- 
ulations by calculating the difference between the true and 
the reconstructed primary energy (applying CIC to the sim- 
ulated data) and was found to be 22% for E =10 17 eV (for 
all primaries) with a slight decrease with increasing energy. 

5 The correction based on a response matrix 

In the reconstruction of an observable, the final result may 
be under- or over-estimating the true value. When repre- 
senting the energy flux as a histogram a particular recon- 
structed event may be stored in the wrong (neighboring) 
energy bin. Thus in every energy bin of our spectrum we 
will have the data correctly belonging to that bin, but also 
data that was migrating from neighboring bins. As the 
energy spectrum is very steep (spectral index 7 w-3) we 
expect that for a given energy bin, the mis-reconstructed 
events falling into it will be coming predominantly from 
lower energy bins. This induces a shift of the reconstructed 
spectrum and a change of the spectral index (due to the 
variable energy resolution of the method with the energy). 
It is possible to account for the effect of fluctuations by cal- 
culating (from simulations) how many events migrate and 
where they add a contribution. Therefore a correction pro- 
cedure based on a response matrix is derived and applied 
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Figure 5: Comparison between the true (known) energy of 
primaries (five masses in fairly equal proportions) used in 
the simulated shower sample and the result of the recon- 
structions: the described reconstruction based on S(500) 
and the result of the standard KASCADE-Grande approach 
based on the N c /j - N M approach; the plot is a profile his- 
togram, the error bars show the spread of data and the con- 
tinuous line is the identity. 
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Figure 6: Reconstructed experimental energy spectrum by 
KASCADE-Grande from S(500)/CIC, multiplied by E 2 5 
together with the result of the standard reconstruction pro- 
cedure (based on the -N M correlation) and the results 
of KASCADE towards lower energies; the continuous lines 
above and below the spectrum show the systematic uncer- 
tainties. 

to the experimental data. The spectrum presented in Fig. 6 
includes the result of this correction. 

6 Conclusions 

The primary energy spectrum has been reconstructed from 
the particular case of charged particle densities recorded 
in the stations of the KASCADE-Grande array at 500 m 
distance from the shower axis, S(500). The study was 
performed for experimentally recorded events and also for 
simulated ones. The analysis of simulated events was 
needed in order to evaluate the reconstruction efficiency 
and quality and to derive a calibration curve Eq - S(500). 
The CIC method was applied on each shower in order 



to correct for attenuation effects. Using the simulation- 
derived calibration the S(500) values are converted into 
primary energy. The S(500)-derived KASCADE-Grande 
spectrum is composition independent. An evaluation of the 
various uncertainty sources has been done and a correc- 
tion based on a response matrix has been employed to ac- 
count for the effects of the fluctuations on the spectral index 
of the reconstructed energy spectrum. The S(500)-derived 
all-particle primary energy spectrum (Fig. 6) shows a shift 
and a slightly different spectral index compared to an inde- 
pendently performed reconstruction approach based on the 
total number of charged particles and the muon number. 
In an event by event comparison, the profile plot in fig. 5 
shows that when reconstructing simulated events with the 
two approaches there is good agreement between the meth- 
ods. Present and future investigations are directed towards 
understanding the origin of the difference between the ex- 
perimental results of the two methods. A possible source 
for this effect could be the shape of the lateral density dis- 
tribution for simulated events that appears to be different 
from the shape of the experimental ones, although there 
could be additional sources for this shift. 
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Abstract: The KASCADE-Grande experiment is an air-shower ground-based observatory designed to study in detail 
the energy spectrum and composition of primary cosmic rays in the region of 10 16 - 10 18 eV. These analyses are based 
on precise measurements of the charged, electron and muon numbers of the cosmic ray air-showers performed through 
different detector systems which come into play simultaneously in KASCADE-Grande during the data acquisition. Due 
to the quality of the data and the number of air-shower observables at disposal through the experiment the collected 
data proves to be also useful to test hadronic interaction models used for air-shower simulations. In this contribution, 
predictions of the QGSJET II-2, SIBYLL 2.1 and EPOS 1.99 hadronic models are confronted with the KASCADE- 
Grande muon data. Besides, the influence of these models on the all-particle energy spectrum derived from the muon size 
is also investigated. 

Keywords: KASCADE-Grande, hadronic models, muons, simulations, energy spectrum 



1 Introduction 

The interpretation of high-energy cosmic ray data relies on 
the extensive air shower (EAS) simulations in which an im- 
portant source of uncertainty is the description of hadronic 
interactions. The main problem is that for the physics of air 
showers the relevant processes lie in the kinematical region 
of small transverse momenta where QCD cannot be applied 
perturbatively and, in most of the cases, where no data 
is available. In this way, phenomenological models and 
parametrizations of accelerator data at low energies must 



be invoked. These models are later extrapolated to high en- 
ergies where they are used as valuable tools to understand 
the EAS observations at cosmic ray observatories [1]. Ac- 
tually, air shower experiments can also help as laboratories 
to test and improve the modern hadronic interaction models 
in energy regions which at the moment are not reachable in 
particle accelerators. This task requires the measurement 
of as many EAS properties as possible with enough preci- 
sion, such as the ones performed with the KASCADE ex- 
periment [2]. In the energy interval of 10 14 — 10 16 , tests of 
early versions of the hadronic interaction models QGSJET, 
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Figure 1: Left: Mean values of the muon correction function shown against the reconstructed muon number N^. Right: 
Fluctuations for the muon number after using the correction function versus the corrected muon number (N'^) assuming a 
mixed primary composition. MC results for different hadronic interaction models are presented. 



EPOS, DPMJET and SIBYLL have been carried out with 
the KASCADE observatory (see [3] and references therein) 
through detailed studies involving the hadron, electron and 
muon contents of air showers. None of the above mod- 
els was able to describe simultaneously all the KASCADE 
EAS data in the energy region around 1 PeV. Now, with 
the enhanced detector KASCADE-Grande [4], the possi- 
bility to extend these tests to the 10 16 — 10 18 eV interval is 
opened. 

First investigations on this subject have been already per- 
formed with KASCADE-Grande. For example, in [5] it 
has been shown that predictions with QGSJET II about the 
muon production height (H^) distributions for EAS with 
zenith angles below 9 < 18° show a discrepancy com- 
pared to the measured data. On the other hand, in [6] 
measurements on the muon lateral distributions of 
EAS performed with KASCADE-Grande were confronted 
with the expected results of QGSJET II-2 and EPOS 1.99 
finding an agreement between experiment and model pre- 
dictions [6]. Some tests regarding the sensitivity of the 
KASCADE-Grande results to the QGSJET II-2 and EPOS 
1.99 hadronic models have also been done, in particular, 
associated with the composition and the all-particle en- 
ergy spectrum of cosmic rays derived from the and the 
charged particle number (N c h) analyses, respectively [6]. 
EPOS 1.99 favors an abundance of light primary particles 
in the data and predicts a higher flux (~ 33%) in compari- 
son with QGSJET II-2 [6]. 

In this work, the hadronic interaction models QGSJET II- 
2 [7], SIBYLL 2.1 [8] and EPOS vl.99 [9] are tested with 
the KASCADE-Grande muon data. The advantage of using 
muons is that in an air shower these particles undergo less 
atmospheric interactions than the electromagnetic compo- 
nent and, in consequence, they reflect directly the physics 
of the first hadronic interactions in EAS. Other tests and 
comparisons employing the p M distributions and the all- 
particle energy spectrum, derived by combining the infor- 
mation from the N c h — Np observables, are performed in 
[10,11]. 



In KASCADE-Grande, measurements of the total muon 
number in EAS {N^, number of muons greater than 
230 MeV) are performed with an array of 192 x 3.2 m 2 
shielded scintillator detectors belonging to the former 
KASCADE experiment [2]. On the other hand, arrival 
times and charged particle densities, employed for estima- 
tions of the EAS arrival direction, N c h content and core 
position, are measured with an enhancement called the 
Grande array, which is composed by 37 x 10 m 2 plastic 
scintillator detectors scattered on a surface of 0.5 km 2 [4]. 

2 Description of data sets 

For the present analysis, all air shower simulations were 
performed with CORSIKA [12] using Fluka [13] to treat 
hadronic interactions in the low energy regime. At high 
energies, the hadronic interaction models QGSJET II-2, 
SIBYLL 2.1 and EPOS vl.99, subjects of this investiga- 
tion, were employed. The response of the detector was 
simulated with a GEANT 3.21 based code. Sets with spec- 
tral indexes 7 = —2.8, —3, —3.2 were produced for several 
primary particle assumptions: H, He, C, Si, Fe and a mixed 
composition scenario (all single primaries in equal abun- 
dances). 

Selection cuts were applied to both experimental and MC 
data. They were chosen according to MC studies to avoid 
as much as possible the influence of systematic uncertain- 
ties in the measurements of the EAS parameters. The se- 
lected data were composed of events with more than 11 
triggered stations in Grande, shower cores inside a central 
area of 1.52 x 10 5 m 2 and arrival directions confined to 
the zenith angle interval of AO = 0° — 40°. These events 
were registered during stable periods of data acquisition 
and passed successfully the standard reconstruction proce- 
dure of KASCADE-Grande [4]. Additionally, only show- 
ers with log Nfj, > 5.1 were considered for this work. Both 
the experimental and simulated data were analyzed and re- 
constructed with the same algorithms. With the above qual- 
ity cuts, the effective time of observation with KASCADE- 
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Fe 


EPOS 1.99 


445 ± 26 


624 ± 31 


636 ± 37 


459 ± 23 


607 ± 30 


624 ± 31 


476 ± 25 


614 ± 30 


604 ± 30 


QGSJET E 


824 ± 33 


832 ± 31 


690 ± 43 


900 ± 40 


833 ± 31 


693 ± 42 


897 ± 100 


825 ± 50 


750 ± 62 


SIBYLL 2.1 


546 ± 44 


657 ± 29 


681 ± 46 


637 ± 39 


672 ± 29 


688 ± 38 


725 ± 44 


681 ± 29 


699 ± 40 



Table 1 : Muon attenuation lengths extracted from Monte Carlo data. The first column represents the hadronic interaction 
model. The corresponding composition scenario and spectral index, 7, of the MC sample under study are specified at the 
upper lines of the table. 



Grande was equivalent to 1424 days. The threshold for full 
efficiency was found at log 10 « 5.4. 



3 Description of the analysis and results 

To start with, all muon data was corrected for systematic 
uncertainties using muon correction functions derived from 
MC simulations for each hadronic interaction model as- 
suming mixed composition and 7 = —3. The functions 
were parametrized with respect to core position, azimuthal 
and zenithal angles, and muon size. In Fig. 1 the mean 
value of the muon correction function for different hadronic 
interaction models is plotted against the uncorrected N^. In 
general, after correction the systematic error on the muon 
number above threshold is found to be almost independent 
of the corrected muon size, N' , and smaller than 6%. Fluc- 
tuations on N' were also investigated. They are shown in 
Fig. 1 for a scenario with mixed composition and 7 = —3. 
Note that although fluctuations exhibit the same tendency 
independently of the hadronic model under consideration, 
they present slight differences in magnitude, as it is the case 
for the correction functions. Therefore, some differences 
are expected when interpreting the same muon experimen- 
tal data with the hadronic interaction models under consid- 
eration. 

In a second step, to test the hadronic interaction models 
with the KASCADE-Grande muon data, predictions on the 
evolution of the muon content with the arrival zenith an- 
gle of the EAS were confronted with observations. The 
task was done comparing the expected and observed val- 
ues of the muon attenuation length, A M . This quantity 
was extracted by applying the Constant Intensity Cut (CIC) 
method to the data as described in reference [14] but using 
a global fit to the attenuation curves, log 10 N'^(9), with the 
known formula 



N' 



<exp[ 



-X o sec(0)/A, 



(1) 



where Xq = 1023 g/cm is the average atmospheric depth 
for vertical showers and N'® is a normalization parameter 
to be determined for each attenuation curve. The results 
for A M are presented in Tables 1 and 2. Discrepancies be- 
tween the experimental values and the simulation results 
can be observed for the studied models. The differences 
do not disappear when modifying the primary composition 
or spectral index. As a consequence, the predicted evolu- 
tion of the muon component with the zenith angle, N'(0) 





A (j (gW) 


Model 




EPOS 1.99 


1851 ± 142 


QGSJET II 


1383 ± 84 


SIBYLL 2.1 


1443 ± 86 



Table 2: Muon attenuation lengths extracted from 
KASCADE-Grande data under the framework of different 
hadronic interaction models. When comparing MC and ex- 
perimental data it should be understood that the same 
correction function was employed. 

(see equation 1) shows also a disagreement with the obser- 
vations. The percentage of deviation for N'(0) between 
experiment and simulations (mixed composition scenario 
with 7 = — 3) inside the frameworks of QGSJET II-2, 
EPOS 1.99 and SIBYLL 2.1 is presented in Fig. 2. For- 
mula 1 was employed to calculate the curves of Fig. 2 us- 
ing the values for A M of Tables 1 and 2 and normalizing 
data in such a way that the predicted N' at 6 — 0° agrees 
with the measured value. Several factors, which are model 
dependent, may come into play in the observed differences: 
from the predicted muon correction function (see graphs in 
Fig. 1), up to the description of the production, evolution 
and fluctuations of the shower, therefore one should be cau- 
tious when extracting conclusions from these differences. 

To finalize this analysis, the KASCADE-Grande energy 
spectra were reconstructed using QGSJET II-2, EPOS 1 .99 
and SIBYLL 2. 1, under the assumption of a mixed compo- 
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Figure 2: Relative deviation of measured data from model 
predictions for the dependence of the muon content with 
the zenith angle in the framework of three hadronic inter- 
action models. 
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Figure 3: Primary energy as a function of the corrected 
muon number estimated from MC simulations for a mixed 
composition scenario with 7 = — 3 using QGSJET II-2, 
EPOS 1.99 and SIBYLL 2.1. 

sition scenario and a spectral index 7 = — 3, and were com- 
pared with each other to test the influence of the hadronic 
interaction model on the KASCADE-Grande muon data. 
Details about the method of reconstruction can be found 
in [14]. One important difference with respect to refer- 
ence [14] is that in the present work, the energy spectra 
were unfolded using the Gold algorithm [15]. As it was 
the case in [14] the primary energy was calculated using 
a calibration function of the form E = a fJi [N' ll (9 re f)} 13 '' , 
where A^(# re /) is the corrected muon number of the EAS 
that is expected at 9 re f according to the CIC method. 
6 re f = 23.3° is a zenith angle of reference. On the other 
hand, a p and (3^ are model dependent parameters, which 
are determined from a fit to the MC data presented in Fig. 
3. Note that at a fixed energy, the predictions of EPOS 1 .99 
for the mean number of muons in the interval 9 = 20° —26° 
give values higher by 14% compared to QGSJET II-2, and 
21% to SIBYLL 2.1, The reconstructed energy spectra are 
shown in Fig. 4 for the full efficiency region. 

As seen in Fig. 4, the all-particle energy spectrum derived 
with QGSJET 11-2 is lower than the spectra calculated with 
SIBYLL 2.1 but larger than the flux estimated with EPOS 
1.99. In particular, at E = 10 s GeV, the above differ- 
ences are of the order of 23% and 36%, respectively. These 
values are comparable to those derived from reference [6] 
(~ 30%) where the charged particle number is used instead 
to reconstruct the energy spectra within the QGSJET II-2 
and the EPOS 1.99 frameworks. In the latter case, how- 
ever, the EPOS 1.99 fluxes are shifted to higher energies in 
comparison with those estimated with QGSJET II-2. 

4 Conclusions 

Three hadronic interaction models: QGSJET II-2, EPOS 
1.99 and SIBYLL 2.1 were tested in this paper by compar- 
ing their predictions for the attenuation of the muon con- 
tent of EAS in the atmosphere with the observations of 
KASCADE-Grande. It was found that the above hadronic 



Figure 4: Energy spectra derived from KASCADE-Grande 
muon data using several hadronic models and assuming a 
mixed composition scenario with 7 = —3. 

interaction models do not describe this aspect of the mea- 
sured muon data. The sensitivity of the all-particle en- 
ergy spectrum derived from the measurements to the 
hadronic interaction models was investigated. The spec- 
trum reconstructed using EPOS 1.99 is lower (~ 36%) in 
comparison with that from QGSJET II-2. The latter being 
smaller (~ 23%) than the spectrum derived with SIBYLL 
2.1 
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Abstract: The muon component of atmospheric air showers is a very important information in astroparticle physics 
because it is related to the primary particle mass and also because it depends on the hadronic interactions happening in 
the shower development. In this paper, we study the muon densities measured by the KASCADE-Grande experiment and 
illustrate its importance in composition studies and testing of hadronic interaction models. The data analysed here was 
measured by the KASCADE-Grande detector and lies in the 10 16 — 10 18 eV energy range. The measured muon density 
is compared to predictions of EPOS 1.99, QGSJet II and Sibyll 2.1 hadronic interaction models. This paper is an update 
of the results presented in [1]. The paper extends its scope by testing two extra hadronic models Sibyll 2.1 and EPOS 
1.99. 
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1 Introduction 

Cosmic rays with energy range between 10 16 and 10 18 eV 
have the potential to reveal interesting astrophysical phe- 
nomena occurring in the Universe. In this energy range the 
transition between galactic to extragalactic predominance 
of the flux might occur. This transition is expected to cause 
an evolution of the mean mass of primary particles. 

The KASCADE-Grande experiment (see figure 1) has 
been set up to measure primary cosmic rays in this en- 
ergy range in order to help in the understanding of these 



questions. The experiment is located at the Karlsruhe In- 
stitue of Technology, Campus North, Germany, where, be- 
sides the existing KASCADE [2] array, two new detector 
set ups (Grande and Piccolo) have been installed. The ex- 
periment is able to sample different components of exten- 
sive air showers (electromagnetic, muonic and hadronic) 
with high accuracy and covering a surface of 0.5 km 2 . For 
an overview of the actual setup of the KASCADE-Grande 
Experiment see ref. [3]. 

Muons are the messengers of the hadronic interactions of 
the particles in the shower and therefore are a powerful 
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Figure 1: Representation of the KASCADE-Grande detec- 
tors. 
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Figure 2: Distribution of the density of muons for three cuts 
in the total number of electrons. 



tool to determine the primary particle mass and to study 
the hadronic interaction models. 

In this article we present studies of the muon component of 
the showers. This paper is an update of a previous publica- 
tion [1] in which the KASCADE-Grande Collaboration has 
shown the muon densities for the first time. In this contri- 
bution, we analyse two more years of data representing an 
substancial increase in the number of events. In this publi- 
cation two extra hadronic interation models are also tested: 
EPOS 1.99 and Sibyll 2.1. 

2 Reconstruction 

As explained in [1] and repeated here, the main parameters 
used in this study are the density of muons and the total 
number of electrons in the shower for which the reconstruc- 
tion accuracy is going to be discussed below. For the recon- 
struction accuracy of the shower geometry see ref. [4]. 

The density of muons is directly measured by the KAS- 
CADE 622 to 2 scintillators. These detectors are shielded 
by 10 cm of lead and 4 cm of iron, corresponding to 20 
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Figure 3: Distribution of the density of muons for three 
distances from the shower axis. 



radiation lengths and a threshold of 230 MeV for vertical 
muons. The error in the measurement of the energy deposit 
was experimentally determined to be smaller than 10% [2]. 

For each shower, the density of muons is calculated as fol- 
lows. The muon stations are grouped in rings of 20 m dis- 
tance from the shower axis. The sum of the signals mea- 
sured by all muon stations inside each ring is divided by 
the effective detection area of the stations. Therefore the 
muon density as a function of the distance from the shower 
axis is measured in a very direct way. No fitting of lateral 
distributions is needed in these calculations. 

The total number of electrons in the shower is reconstructed 
in a combined way using KASCADE and KASCADE- 
Grande stations. A lateral distribution function (LDF) of 
the Lagutin type can be fitted to the density of muons mea- 
sured by the KASCADE detector [5]. After that, using 
the fitted function, the number of muons at any distance 
from the shower axis can be estimated. The KASCADE- 
Grande stations measure the number of charged particles. 
The number of electrons at each KASCADE-Grande sta- 
tions is determined by subtracting from the measured num- 
ber of charged particles the number of muons estimated 
with the LDF fitted to the KASCADE stations. 

At this stage, the number of electrons at each KASCADE- 
Grande station is known. Finally, a modified NKG [6] 
function is fitted to this data and the total number of elec- 
trons is determined in the fit. 

Quality cuts have been applied to the events in this analysis 
procedure. We have required more than 19 KASCADE- 
Grande stations with signal. The showers used in all anal- 
yses along this paper were reconstructed with zenith angle 
between and 42 degrees. The same quality cuts were ap- 
plied to the simulated events used for reconstruction studies 
and to the data presented in the following section. After the 
quality cuts, the total number of electrons can be estimated 
with a systematic shift smaller than 10% and a statistical 
uncertainty smaller than 20% along the entire range con- 
sidered in this paper [4]. 

Figure 2 shows the measured density of muons at three dis- 
tances from the shower axis for events with a total number 
of electrons (iV e ) in the range 7.0 < Logl0(N e ) < 7.3 
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Figure 4: Measured distribution of the density of muons at Figure 6: Measured distribution of the density of muons at 
400 m compared to the predictions of QGSJet II. 400 m compared to the predictions of EPOS 1 .99. 
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Figure 5: Measured distribution of the density of muons at 
400 m compared to the predictions of Sibyll 2. 1 . 

(« 10 17 eV). Similar plots were obtained for other N e 
ranges. 

Figure 3 shows the density of muons at 400 m from the 
shower axis for events with total number of electrons 
(N e ) in the range 6.7 < Logl0(N e ) < 7.0, 7.0 < 
Logl0(N E ) < 7.3 and 7.3 < Logl0(N e ) < 7.6. Simi- 
lar plots were obtained for other distances from the shower 
axis. 

Figure 2 and Figure 3 show the general expected trend: a) 
decrease of the muon density with increasing distance from 
the shower axis and b) increase of the muon density with 
increasing total number of electrons. In the next sections 
we explore these relations in order to show the capabilities 
of the KASCADE-Grande experiment for a composition 
study and for tests of the hadronic interaction models. 

We present data for 7.0 < Logl0(N e ) < 7.3 and 390 < 
r < 410 to, these cuts have been chosen in order to min- 
imize the fluctuation of the signal and the reconstruction 
inaccuracy and to maximize the number of showers for 
which we have data, however the same conclusions would 
be drawn for all parameter cuts. 

3 Simulation 

For all studies in this paper we have used the CORSIKA 
[7] simulation program with the FLUKA [8] option for low 
energy hadronic interactions. Three high energy hadronic 



interaction models were used EPOS 1.99 [9], QGSJet II 
[10] and Sibyll 2.1 [11]. No thinning is used [7]. 

CORSIKA showers are simulated through the detectors and 
reconstructed in the same way as the measured data, such 
that a direct comparison between data and simulation is 
possible. 

Figures 4, 6 and 5 show the comparison of the measured 
density of muons to values predicted by QGSJet II, EPOS 
1.99 and Sibyll 2.1, respectively. For the three hadronic in- 
teractions models we show the limiting cases of proton and 
iron nuclei as primary particles. It can be seen in figures 4, 
6 and 5 that the data lie well within the proton and iron lim- 
its for QGSJet II, EPOS 1 .99 and Sibyll 2.1. These graphics 
are going to be further discussed in the next sections. 

4 Analysis 

Figure 7 shows the mean muon density as a function of the 
distance from the shower axis compared to the predictions 
of QGSJet II, EPOS 1 .99 and Sibyll 2. 1 . The three hadronic 
interaction models bracket the data within the proton and 
iron limits for the entire range of distances from 100 to 720 
meters. For distances further than 720 meters the statistics 
is not enough for a conclusion. 

Interesting to note is also the slope of the LDF Consider- 
ing an equal probability trigger for protons and iron pri- 
maries as a function of distance from the shower axis, one 
can analyse the slope of the LDF as a test of the hadronic 
interaction models or as a composition parameters. 

The intermediate distance from 200 to 600 m is the most 
significant range. For larger distance the statistical fluctu- 
ations dominates and at short distances the accuracy of the 
core position reconstruction impacts the measurements of 
the LDF resulting in a steeply falling curve as seen in figure 
7. 

It is clear that the measured LDF is not parallel to the 
QGSJet, nor EPOS 1.99 nor Sibyll 2.1 curves. That shows 
that the slope of the LDF cannot be well described by nei- 
ther models if the composition is supposed to be dominated 
by only one light or only one heavy primary. If the compo- 
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Figure 7: Lateral distribution of muons compared to the 
predictions of QGSJet II, EPOS 1.99 and Sibyll 2.1. 

sition is mixed the slope varies with distance and the inter- 
pretation depends on the exact composition. 

Figure 8 shows the evolution of the mean muon density 
as a function of N e , The calculations done with QGSJet 
II, EPOS 1.99 and Sibyll 2.1 using proton and iron nuclei 
as primary particles bracket the data in the entire range of 

5 < Logl0(N e ) < 8. 

Nevertheless, both figures 7 and 8 show that EPOS 1.99 
would require a very light primary composition in order to 
fit the data. On the other hand, QGSJet II and Sibyll 2.1 
could fit the data with an intermediate primary abundance 
between proton and iron nuclei. 

Besides that, in figure 8 it is possible to analyse a possible 
transition of the primary component with increasing total 
number of electrons. The analysis done with both mod- 
els show no abrupt change in the compositon in the entire 
energy range. 

5 Conclusions 

The Grande array is in continuous and stable data taking 
since December 2003. The quality of the detector can be 
illustrated by the smooth data curve and small fluctuations 
in figures 7 and 8. 

In this article, we have briefly described the procedure used 
to measure the density of muons with the KASCADE array 
and we have studied its correlation with the distance from 
the shower axis and the total number of electrons in the 
shower. 

The density of muons in the shower is measured directly by 
the KASCADE detectors. We have used this data to study 
the hadronic interaction models QGSJet II, EPOS 1.99 and 
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Figure 8: Muon density as a function of the total number of 
electrons compared to the predictions of QGSJet II, EPOS 
1.99 and Sibyll 2.1. 

Sibyll 2.1. The plots shown here are an update of the pre- 
vious publication [1] 

Figure 8 shows no abrupt change with increasing total num- 
ber of electrons up to LoglQ(N e ) = 7.5 « 5 x 10 17 eV. 
The mean primary mass estimation would depend on the 
hadronic interaction model used. 
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Abstract: With the Muon Tracking Detector in the KASCADE-Grande experiment mean EAS muon pseudorapidities 
are investigated. Here we report on the results of studying the sensitivity of this quantity to the mass of primary cosmic 
ray particles. Obtained values of the mean logarithmic mass in the 10 16 eV - 10 17 eV range of primary energies, based on 
the QGSIetll - FLUKA interaction model combination, are compared with the results of other experiments. The validity 
of the model in reproducing experimentally measured pseudorapidity values and its comparison with the EPOS 1.99 is 
discussed. 
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1 Introduction 

The Muon Tracking Detector (MTD) [1] is one of the de- 
tector components in the KASCADE-Grande EAS exper- 
iment [2] (see Fig.l), operated at the Karlsruhe Institute 
of Technology (KIT) - Campus North, in Germany, by an 
international collaboration. The MTD measures directions 
of muon tracks in EAS with excellent angular resolution 
of w 0.35°. These directional data allow to investigate 
the longitudinal development of the muonic component in 
showers which is a signature of the development of the 
hadronic EAS core, being in turn dependent on the mass of 
the primary cosmic ray particle initiating a shower. Such 



studies can be done either by the determination of a mean 
muon production height [3] or by using the mean pseudo- 
rapidity (r/) of EAS muons, expressed in terms of their tan- 
gential (r) and radial (p) angles (quantities reconstructed 
in the experiment) [4], [5]. In this work we investigate to 
what extent one can use the muon pseudorapidity for the 
determination of primary mass. 

2 Muon pseudorapidities in KASCADE- 
Grande 

In KASCADE-Grande muons can be registered up to 
700 m from the shower core, but normally different anal- 
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Figure 1 : Layout of the KASCADE-Grande experiment distributed over the KIT - Campus North area. KASCADE is 
situated in the North-East corner of the Campus; note the position of the Muon Tracking Detector. 



yses are carried out in specific distance ranges. In Fig. 
2 lateral distribution of mean EAS muon pseudorapidities 
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Figure 2: Lateral distribution of mean muon pseudora- 
pidity in the limited distance range (see text) measured 
in KASCADE-Grande and compared with the CORSIKA 
simulation values for proton and iron primaries. 

measured with the MTD in KASCADE-Grande are shown, 
together with the results of CORSIKA [6] Monte-Carlo 
simulations, employing a QGSJetll [7] and FLUKA2006 
[8] model combination. Requirement of shower size 
logN e > 6 ensures full experiment trigger efficiency 
for registration, both, proton- and iron-initiated show- 
ers. Muons from iron showers are, per average, produced 
higher than those from proton showers, and this leads to 
the observed difference in the mean pseudorapidities on 



ground, i.e., to the primary mass sensitivity of muon pseu- 
dorapidities. In Fig. 2 one can notice in the experimental 
data, below 200 m, a bias towards proton mass, resulting 
(for showers above 10 16 eV) from some saturation effects 
in the MTD close to the shower core. Above 300 m, an 
increasing bias towards iron mass is observed, caused by 
non-equal MTD trigger efficiency for all types of primaries 
there. Therefore, all subsequent investigations are done in 
the distance range 200 m - 300 m, where, as it was checked 
with the N^/Ne ratio of showers used in the analysis, regis- 
tration efficiency of the MTD is primary mass independent. 

3 Mean logarithmic mass calculated with 
muon pseudorapidities 

The distribution of measured muon pseudorapidity in the 
range 200-300 m from the shower core shows, that its mean 
value is contained between the simulated values for proton 
and iron primaries (Fig. 3). 

In order to assess the average mass composition of cosmic 
rays with the standard < InA > procedure one has to be 
sure, that the mean muon pseudorapidity linearly depends 
on primary mass. The results of the linearity check are 
shown in Fig. 4. Here, the < InA > was calculated from 
n distribution for simulated carbon showers and compared 
with the known value lnl2=2.49. The mean value of calcu- 
lated InA differs here by less than 2% from the true value 
for carbon primary, justifying the use of muon pseudora- 
pidity for the determination of < InA > of cosmic rays 
above 10 16 eV. 
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log [Eo/GeV] 


< E > [GeV] 


< InA > 


< InA > 


range 




(QGSJetll-FLUKA) 


EPOS 1 .99-FLUKA 


7.0 - 7.3 


(1.38±0.01)xl0 7 


0.4± 0.3 


1.0± 0.5 


> 7.0 


(2.23±0.01)xl0 7 


1.1± 0.2 


1.2± 0.3 


7.3 - 7.6 


(2.69±0.01)xl0 7 


1.6±0.3 




> 7.3 


(3.92±0.03)xl0 7 


2.0± 0.2 


1.5±0.5 


7.6 - 7.9 


(5.34±0.02)xl0 7 


2.2±0.5 


2.4± 1.3 


> 7.6 


(7.12±0.08)xl0 7 


2.7±0.4 


3.5± 1.1 


> 7.9 


(13.20±0.26)xl0 7 


2.5± 1.4 





Table 1: Results of the < InA > calculated in the 10 16 — 10 17 eV primary energy range for two high-energy interaction 
models. 
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Figure 3: Distribution of measured muon pseudorapidity 
in the range 200-300 m from the shower core together with 
the simulation results for proton and iron initiated showers. 
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Figure 4: Test of the linear dependence of mean EAS muon 
pseudorapidity on the logarithm of primary mass. Triangles 
- simulated < n > values for H,C and Fe versus true InA on 
the x-axis. A full circle - < InA > calculated for carbon, 
using its < n > value. 



In the following the results for the primary energy decade 
10 16 — 10 17 eV will be given. The primary energy was 
calculated using total number of charged particles and total 
number of muons in showers, with a procedure described 
in [9]. The energy ranges in which the < InA > was calcu- 
lated are shown in the first column of Table 1 . In the second 
column mean energy values, taking into account the spec- 
trum with an index -3.1, are given. Columns 3 and 4 con- 
tain calculated < InA > values for QGSJetll and EPOS 
1.99 [10], respectively. Due to the very limited statistics of 
available simulations with EPOS not in all energy bins the 
calculation was possible. In those, where it was possible, 
the errors (only statistical were considered) are, anyway, 
much larger than in case of the QGSJetll. 

The results from columns 3 and 4 of Table 1 are presented 
in Fig. 5, together with a collection of < InA > values 
(taken from Fig. 14 (top) in Ref. [11]) derived from the 
average depth of the shower maximum by various experi- 
ments. 



4 Discussion and conclusions 

As seen in Fig. 5 our results, obtained with EAS muons 
pseudorapidity measured on ground level, are compatible 
with the collection of mean logarithmic mass values de- 
rived from the average depth of the shower maximum. 
They show similar rise of < InA > with energy in the in- 
vestigated primary energy range. This compatibility, rather 
with the results obtained from measurement of the light 
generated during the shower development than the ones 
measured on ground (e.g. electron/muon ratio) is because 
the mean muon pseudorapidity is also a signature of shower 
development. 

Obtained values are generally lower than those from other 
experiments. The main reason for this can be, that the 
QGSJetll model, despite improvement with respect to 
QGSJetOl, gives still too high mean pseudorapidity values 
compared with the measurements. 

It was found in previous studies (e.g. [12]), that in the in- 
vestigated distance range from the core about 40 % of the 
registered muons is produced in the hadronic interactions 
above 200 GeV, i.e., they are modeled in the simulations 
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Figure 5: Results of mean logarithmic mass values obtained with mean muon pseudorapidities, shown together with the 
values derived from the average depth of the shower maximum by different experiments, from Fig. 14 (top) of Ref. [11]. 



by QGSJetll. And this model creates more muons at the 
heights above 4-5 km, than observed experimentally [3]. 
Such muons have higher mean muon pseudorapidity than 
the ones created deeper in the atmosphere [13], thus affect- 
ing the determination of the < In A >. 

The mean pseudorapidity values obtained in simulations 
suffer from the insufficient shower statistics in the simu- 
lations. A standard approach is to use the same simulated 
shower 5 or 10 times over the detector area in order to re- 
duce the time needed for simulations to acceptable value. 
In this particular investigation such an approach creates ad- 
ditional, systematic error in the determination of < In A >, 
impossible to be quantified without using each simulated 
shower once, what would require very significant increase 
of the simulation time. 

In Fig. 5 some < InA > values obtained with the EPOS 
1 .99 model are shown for the comparison. Due to the much 
smaller number of simulated data the statistical errors are 
much larger than in case of QGSJetll (see also Table 1). 
Moreover, the above mentioned effect of shower statistics 
is influencing the results even more. 

One can say that the general behaviour showing rise of the 
primary mass with the energy is also seen with this model, 
however it is not possible to compare the values of the log- 
arithmic mass predicted by the two models in question. 

In conclusion, the mean muon pseudorapidity in EAS is a 
primary mass sensitive parameter and, providing the large 
enough number of simulated showers is available for the 
analysis, the < InA > of primay cosmic rays can be de- 
rived. However, the small difference between < n > values 
for proton and iron initiated showers (w 0.2) would require 
a significant increase in the number of simulated showers 
and, even then, the remaining uncertainties would prohibit 



to make a convincing mass composition studies with, at 
least, 3 mass groups. 
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Abstract: A large area (128m 2 ) streamer tube detector, located within the KASCADE-Grande experiment, has been built 
with the aim to identify muons and their directions from extensive air showers by track measurements. We discuss the 
possibility of observation of Gamma-Ray sources by means of single isolated muons above the background of cosmic- 
ray muons using a muon tracking detector (MTD) exhibiting good angular resolution. Properties of the pion photo- 
production process and of the MTD which support the identification of Gammas are discussed. Preliminary Gamma 
spectrum accumulated from Crab and the Mkn421 flux correlation with X-ray (RXTE/PCA) are presented. 

Keywords: gamma ray sources, muon tracking 



1 Introduction 

A reliable understanding of the muon production by pri- 
mary gamma-rays is mandatory to gauge the sensitivity of 
the experiment to gamma primaries. The Vector Meson 
Dominance Model is usually employed for photo-nuclear 
interactions at gamma energies above a few GeV To il- 
lustrate the critical role played by the event generators in 
predicting the muon content of showers, the Feynman-x 
distribution of charged pions (in the laboratory frame) as 
calculated by the FLUKA code demonstrates [1] a basic 
difference between gamma and proton-induced collisions: 
The gamma primaries lead to a much larger fraction of 



high-x secondaries than the proton primaries, therefore, fo- 
cussing the pions to very forward direction, and, therefore, 
compensating for the much smaller pion photo-production 
cross section. 

High energy gamma rays produce muons in the Earth's at- 
mosphere that can be detected and reconstructed in relative 
shallow underground muon detectors. Such detectors are 
sensitive to muon energies of a few GeV. Although muons 
of such low energy compete with a large background of 
cosmic ray muons, they can be identified provided the de- 
tector has sufficient effective area and resolution. Unlike 
air-Cherenkov telescopes [2, 3, 4] muon detectors cover a 
large fraction of the sky with a large duty cycle. The advan- 
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tage is considerable in studying the emission from highly 
variable sources. Moreover, background multi-muon bun- 
dles can be conveniently rejected without suppression of 
the predominantly single-muon gamma signal. 

Detected muons originate in gamma induced pion produc- 
tion with some 20 times higher gamma energies [5]. A 
low energy threshold for muon detection provides for the 
correspondingly low energy gammas a deep view into the 
Universe. Muons of > 1 GeV associated with primary 
photons of several GeV provide an almost attenuation-free 
window into the depth of the universe suffering little from 
flux losses due to collisions with IR and CMB radiation via 
pair production. 

2 Muons from Gammas 

Gamma rays initiate atmospheric cascades of mostly elec- 
trons and photons, but also some muons. Muons originate 
from the decay of charged pions which are also photo- 
produced by high energy shower photons [6] albeit with 
much smaller cross section. 

The signal-to-noise ratio, defined as the number of events 
divided by the square root of the number of background 
events in a resolution pixel of a° x c , depends on the de- 
tector area A and the zenith angle as [6]: S/N 1 / 2 = 
A 1 / 2 /cos@ 9 <7. The formula simply expresses that the 
signal-to-noise ratio is improved for increased area A, bet- 
ter resolution a and sources observed at larger zenith angle 
0, where the cosmic ray background muon rate is reduced. 

The MTD [7] is sensitive to the gamma energy region 
above ~ 10 GeV while the muon energy cut equals to 
0.8 GeV. 

As mentioned above the background includes some frac- 
tion of multi-muon events. Rejecting multi-muon events 
not only improves the signal-to-noise ratio, it also improves 
the angular resolution which may be degraded by less reli- 
able reconstruction of complex muon bundles initiated by 
high energy cosmic ray particles which are accompanied 
by other shower particles. MTD with its detection area of 
Amtd = 128m 2 has below the shielding an average rate 
of 2.5kHz which results above an energy of 0.8 GeV in 
~ 10 7 tracks from background muons per year in a 1° x 1° 
pixel. Assuming a further reduction because of the strong 
focusing of the gamma induced muons by about 10, this 
rate may lead for Crab [6] to S/N 1 ' 2 = 40 for one year of 
running. 

3 Muon Tracking Detector (MTD) 

The MTD [7] is built out of Streamer Tube (ST) chambers 
of 4 m length and located in the KASCADE-Grande experi- 
ment [9]. The ST chambers are grouped in, so called, mod- 
ules. Four modules, three positioned on horizontal planes 
(top, middle, bottom) and one arranged vertically (wall), 
form a muon telescope. The whole detector comprises 16 



telescopes arranged in two rows. The low mass structure 
of the detector design reduces secondary interactions in the 
sensitive detector part. 

To deal with clean muon tracks an efficient filter ab- 
sorbing a large fraction of the low energy electromag- 
netic component of proton induced showers is manda- 
tory. The larger Bethe-Heitler cross-section induced show- 
ers for gammas leading to pair production and subsequent 
electromagnetic cascade are readily absorbed by the fil- 
ter. The energy threshold (0.8 GeV) is suitable for back- 
ground reduction. Following the characteristic energy e M = 
m^QAkm/r^c ~ 0.8 GeV muons from proton induced 
showers do not survive over the path length of 6.4 km and 
with energies close to the MTD threshold but decay to elec- 
trons which are readily absorbed. Photo-nuclear produced 
muons survive more easily because of comparatively lower 
production height. This situation improves the S/N ratio 
for photon detection. The stability of the MTD [7] is very 
good. The detector gas system follows precisely the atmo- 
spheric gas pressure. 

The first-level 7/proton discrimination is mostly achieved 
by the characteristics of the MTD. The ratio of the number 
of hits in the top module to the track number Nf°f / N tra ck 
in each tower is a powerful tool in discriminating against 
hadron induced air showers. Only tracks with one hit in 
each module are accepted. Only hits with small cluster 
size 'els' (clusters of readout wires or readout strips) are 
accepted. The background of high energy electrons or 
hadrons is further reduced by cuts in the track quality Q 2 
(see [7]). Showering electrons or hadrons lead to larger 
'els' and, therefore, smaller Q 2 . 

For the subsequent analysis the following cuts in the data 
have to be employed: 1) Only data over one full sidereal 
day (86163 solar seconds) are considered. 2) Only sidereal 
days where all 16 telescopes are functioning are used. 3) 
The daily rate should follow a Gaussian distribution, and 
only days with rate within ±4cr are accepted. 4) Only data 
corrected for pressure and temperature in the atmosphere 
and in the detector are included. 

4 Gamma- Source Search 

To identify sources which emit high energy gammas and 
which are with small probability converted to single muons 
deep in the atmosphere, we have the possibility to restrict 
the arrival direction and arrival time. The restriction in 
time and direction has the task to reduce other disturbing 
sources like the Sun. We have to cope with about 3 orders 
of magnitude larger background from cosmic ray muons in 
a canonical 1° x 1° window in the sky. 

The optimal square bin size for the search of a point 
gamma-ray source with the MTD is 1° on a side, corre- 
sponding to a Gaussian angular resolution of about 0.3°. 

For the search of point sources, the fluctuation of events in a 
fixed direction of the sky is investigated. When examining 
the fluctuation of the number of muon events bin by bin and 
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searching for possible sources with a muon event excess, 
we employ the numbers N on , T on , N ff, T a ff. We use Li- 
Ma [10] formula to examine the fluctuation of the number 
of muon events bin by bin and to search for possible muon 
event excess. In the current analysis this procedure is used 
with bin size AS = 0.25° and Aa = AS/cosS, confined to 
specific regions on the sky to investigate the MTD sensitiv- 
ity to specific gamma-ray sources. Only specific track di- 
rections in and $ and only specific arrival time intervals, 
where the corresponding source is high (winter months for 
Crab) (Q source < 35°) are chosen to accumulate clean sin- 
gle muon events in the S versus a (declination versus right 
ascension) plane. 

The analysis shows that the strong requirements of only 
1 hit in each detector module and all els < 3 give good 
profiles for Crab and Mkn42 1 . 

Observation of Mkn421 with the ARGO-YBJ experiment 
is reported by S.Vernetto et al. [11]. The FERMI view of 
the TeV blazar Mkn421 is given in [12] and the MAGIC 
observations of Mkn421 and related optical/X-ray /TeV 
multiwavelength studies are reported in [13]. Mkn421 
is a very active Galaxy nucleus showing strong variation 
in the X-ray and TeV Gamma-ray fluxes and appears at 
38.3°, 11.08 h which is 0.1° closer to the zenith and 0.15° 
later in time than the nominal position. The muon inten- 
sity is sampled and weighted with the actual MTD effi- 
ciency and smoothed afterwards. Again, only data epochs 
are considered when Mkn421 is high. Interestingly, the 
source profile is improved when employing days (prop- 
erly matched in time) for which appreciable X-ray flux 
(> 2.5 photons / day) is reported [15]. For the Mkn421 
source a correlation analysis with the X-ray flux measure- 
ments by the All-Sky-Monitor (ASM) onboard the RXTE- 
satellite [15] was performed providing the correlation as 
shown in Fig. 1. The picture exhibits a clear correlation 
between the time averaged (1 day) high energy muon flux 
and the time averaged X-ray flux above 2 keV. A narrow 
window around the nominal Mkn421 position improves the 
correlation in Fig. 1. Also a correlation between the muon 
flux from the MTD in the region of Mkn42 1 and the gamma 
flux in the GeV range recorded by FERMI satellite [16] 
is observed. The dashed line in Fig. 1 represents a fit to 
the correlation and deserves further studies. Employing 
days with almost constant X-ray activity from Mkn421, the 
conversion of GeV gammas into muons in the atmosphere 
above the MTD can be investigated. Fig. 2 shows the cor- 
relation of the muon rate with the atmospheric pressure 
overburden. It shows the strong gain of photo-produced 
pions with pressure especially after correction (solid points 
in Fig. 2) for the pressure dependent track detection effi- 
cieny of the MTD [7]. The strong dependence of the muon 
rate on the atmospheric pressure suggests the pion photo- 
production to occur in the lower atmosphere. 

Fig. 3 shows for Crab a two-dimensional multi-quadric in- 
terpolation [8] of the data points in the <5 versus a (decli- 
nation versus right ascension) plane (in the boundaries of 
a MAGIC presentation). The contour scale is in units of 
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Figure 1 : Correlation of the X-ray flux with the muon rates 
in a narrow window around the nominal Mkn421 position. 
Data were taken in a high flux period during February 2010. 
The dashed line emphasizes a linear dependence between 
the X-ray flux and the GeV gamma ray flux. So far, only 
fraction of all data are considered. 

excess counts per smoothing radius. The position of the 
source as quoted by the MAGIC group [14] is given by 
the cross. The displacement of the source to larger zenith 
9, and smaller azimuth $, and possible deformation is due 
to the geomagnetic field B geo and is expected to depend 
on the momentum of the muons and their orientation with 
respect to the geomagnetic field components. 

Gamma-ray induced muons in the deep atmosphere we ex- 
pect to be almost charge symmetric. Therefore, the spread 
of the source profile is reproducing the muon momentum 
distribution. Muons are considered to exhibit negative 
charge if they are deflected towards the East from the line 
connecting the actual Crab position with the geomagnetic 
field direction. Mean deflection of about 0.6° is expected 
for 1 GeV muons created 1 km above ground. Only nega- 
tively charged muons to the East side of the source are con- 
sidered for accumulating an angle-distance spectrum. The 
frequency of muon tracks as function of 1 / distance to the 
nominal Crab source position is plotted in the Fig. 4. result- 
ing in steep falling spectrum. The inverse angle-distance is 
expressed in muon momentum, assuming a ~ p^ x B geo 
dependence. Assuming further that muons at threshold 
(0.8 GeV) stem from gammas about 10 times higher in mo- 
mentum would provide a preliminary gamma energy scale. 
The preliminary flux normalization considers the gamma 
flux attenuation in the air, the small tr 7X cross section on 
air, the mass density of the atmosphere and the momentum 
dependence of the n =$> /i conversion rate and the exposure 
of the MTD. The dashed line represents a flux^ ~ p~ 2 
flux dependence. A more detailed unfolding of the gamma 
source spectrum has still to be carried out. 
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Figure 2: Variation of the photo-pion production yielding 
muons with atmospheric pressure for a window of almost 
constant X-ray flux from Mkn421 (see Fig. 1, dashed re- 
gion). The full symbols consider the effect of the MTD 
track detection efficiency. 
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Figure 3: Muon yield close to the nominal Crab position 
in a 3° x 3° field of view. The cross represents the Crab 
position after MAGIC [14]. The contour lines vary in 25% 
steps. So far, only fraction of all data are considered. 
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5 Outlook 

The high resolution muon tracking detector MTD in the 
KASCADE-Grande experiment demonstrates to be capa- 
ble to identify gamma point sources in the sensitivity range 
of Crab fluxes. The highly variable gamma source Mkn421 
provides in its 'high' state a test beam to study the effi- 
ciency of 7 =>■ fj, conversion in the atmosphere depend- 
ing on the atmospheric parameters, and further tuning of 
the MTD response to gammas. Future analysis of a larger 
data sample will provide more detailed information on the 
nature of high energy gamma source muons. There is a 
common understanding that the high energy gamma source 
muons serve as sensitive probes to investigate the high en- 
ergy photon interactions in the atmosphere, providing a 
gamma detector in the multi-GeV range. 
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Figure 4: Muon flux variation derived as function of the 
inverse distance from the nominal source position and ex- 
pressed in momentum units above the muon momentum 
threshold of the MTD (0.8 GeV). So far, only fraction of 
all data are considered. 



